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1.  Objective : @3 6 . 1@UNARR. 


OB  JT . 


I  !The  objective  of  this  research  was  to  explore  new  concepts  for 


boundary  layer  control  using  neural  networks  and  surface  arrays  of  sensors 


and  actuators . 


2  .  Approach :  @3  7 . 1@UNARR .  -  APPR . 

The  research  conducted  numerical  simulations  of  turbulent  boundary  layers 
under  the  influence  of  distributed  arrays  of  nerual  network  controlled  sensors 
and  actuators . 


3.  Progress : @3 8 .1@UNARR.  -  PROG.  -  FROM  DD  MMM  YY  TO  DD  MMM  YY 

Insertion  of  a  microflap  into  the  flow  domain  creates  an  internal  boundary 
which  must  be  simulated  properly.  We  used  an  idea  of  Peskin  and  McQueen  (1989) 
of  using  a  body  force  to  represent  the  dynamical  effect  on  the  fluid  of  a 
moving  wall  immersed  in  the  fluid.  In  their  method,  and  our  applicationl  the 
position  of  the  internal  boundary  within  the  flow  domain  does  not  generally 
coincide  with  the  computational  mesh.  Flow  properties  are  interpolated  from 


the  computational  mesh 

to 

the  locations 

of 

boundary 

nodes  in 

order 

to 

calculate  the  force 

of 

the 

boundary  on 

the 

fluid. 

Then  that 

force 

is 

distributed  onto  nodes 

of 

the 

computational 

mesh 

near  the 

internal  boundary . 

The  force  in  Peskin  and  McQueen's  work  was  based  on  a  balance  of  tension  in  an 


elastic  wall  against  the  hydrodynamic  forces.  We  instead  use  an  idea  developed 
by  Jamaludin  Mohd-Yusof  (1997)  of  creating  a  fictious  force  constructed  so  as 
to  drive  the  velocity  at  desired  locations  to  desired  values.  In  his 
application,  the  desired  locations  are  a  riblet  surface  and  the  desired  values 
are  zero  to  fulfill  the  no-slip  condition.  In  our  application,  the  desired 
locations  are  the  flap  positions  and  the  desired  values  are  the  flap 
velocities,  also  to  fulfill  the  no-slip  condition.  A  difference  in  Mohd- 
Yusof'  s  application  and  ours  is  the  moving  boundary  created  by  the  flap. 
Another  difference  is  that  the  flap  is  an  internal  boundary  with  flow  on  both 
sides.  Mohd-Yusof' s  riblet  surface  is  an  external  boundary  with  only  one  side 
within  the  flow  domain.  There  was  some  question  whether  we  could  use  this 
method  with  an  internal  boundary.  It  was  used,  however,  to  simulate  internal 
boundaries  in  some  simple  flows  with  known  analytic  solutions,  and  it  did 
converge  to  the  proper  results. 

4 .  Progress : 

A  series  of  numerical  experiments  were  performed  to  explore  new  control 
strategies  suitable  for  distributed  control.  All  computations  were  conducted 
in  two-dimensional  channel  geometry.  The  Reynolds  numbers,  based  on  the 
wall-shear  velocity  and  channel  half-height,  were  between  100  and  400.  A 
spectral  channel  code  ,was  used  for  all  computations.  A  summary  of  various 
control  schemes  explored  is  given  along  with  relevant  publications,  to  which 
the  interested  reader  is  referred  for  detailed  descriptions  of  each  control 
algorithm.  One  of  these  control  algorithms  (a  simple  feedback  law  obtained 
from  a  neural  network)  was  implemented  into  an  M3  chip.  Detailed  descriptions 
of  the  M3  chip  will  be  reported  in  the  final  report  of  the  AFOSR/DARPA  program 


mentioned  above. 
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1.  Background  and  Objective 


This  research  effort  was  conducted  in  close  collaboration  with  an  AFOSR/URI  grant  (F49620-93-1- 
0332,  Integrated  Microtransducer  and  Neural  Network  Systems  for  Distributed  Control),  an 
AFOSR/DARPA  grant  (F49620-97-1-0514,  Distributed  Turbulent  Flow  Control  by  MEMS 
Integrated  with  Neural  Networks),  and  an  AFOSR  grant  (F49620-95- 1-0263,  Advanced  Concepts 
for  Boundary  Layer  Control).  The  objectives  of  the  first  two  programs  were  to  develop  a 
microsensor-microelectronics-microactuator  (M3)  system  on  a  chip,  which  can  be  used  in  real-time 
control  for  reducing  the  viscous  drag  in  turbulent  boundary  layers.  The  third  AFOSR  grant 
mentioned  above  and  the  present  grant  were  awarded  to  assist  the  design  of  microelectronics  by 
developing  new  control  strategies  suitable  for  the  M3  system.  Since  the  control  logic  must  be 
included  on  a  chip  that  contains  multiple  sensors  and  actuators,  it  was  imperative  to  develop  a 
simple  control  algorithm.  Various  different  approaches,  ranging  from  a  nonlinear  adaptive 
controller  based  on  a  neural  network  to  a  linear  controller  based  on  a  linear  systems  theory,  were 
explored  by  testing  in  numerical  experiments. 

2.  Accomplishments 

A  series  of  numerical  experiments  were  performed  to  explore  new  control  strategies  suitable  for 
distributed  control.  All  computations  were  conducted  in  two-dimensional  channel  geometry.  The 
Reynolds  numbers,  based  on  the  wall-shear  velocity  and  channel  half-height,  were  between  100 
and  400.  A  spectral  channel  code  was  used  for  all  computations.  A  summary  of  various  control 
schemes  explored  is  given  below  along  with  relevant  publications,  to  which  the  interested  reader  is 
referred  for  detailed  descriptions  of  each  control  algorithm.  One  of  these  control  algorithms  (a 
simple  feedback  law  obtained  from  a  neural  network)  was  implemented  into  an  M3  chip.  Detailed 
descriptions  of  the  M3  chip  will  be  reported  in  the  final  report  of  the  AFOSR/DARPA  program 
mentioned  above. 


2.1  Neural  Network 

The  first  control  scheme  considered  was  an  adaptive  controller  based  on  a  neural  network. 
Motivation  of  this  work  was  a  successful  demonstration  of  drag  reduction  from  the  opposition 
control  by  Choi  et  al.  (1993).  The  major  drawback  of  the  opposition  control,  however,  is  that  it 
requires  information  inside  the  flow,  rendering  the  approach  impractical.  A  neural  network  was 
constructed  and  trained  as  an  inverse  model  of  a  turbulent  channel  flow.  The  inverse  model  was 
used  as  an  adaptive  controller.  A  simple  control  network,  which  employs  blowing  and  suction  at 
the  wall  based  only  on  the  wall-shear  stresses  in  the  spanwise  direction,  was  shown  to  reduce  the 
skin  friction  by  as  much  as  20%  in  direct  numerical  simulations  of  a  turbulent  channel  flow.  Also,  a 
stable  pattern  was  observed  in  the  distribution  of  weights  associated  with  the  neural  network.  This 
allowed  us  to  derive  a  simple  control  scheme  that  produced  the  same  amount  of  drag  reduction. 
This  simple  control  scheme  generates  optimum  wall  blowing  and  suction  proportional  to  a  local 
sum  of  the  wall-shear  stress  in  the  spanwise  direction.  The  distribution  of  corresponding  weights  is 
simple  and  localized,  thus  making  real  implementation  relatively  easy.  Further  details  on  this  work 
can  be  found  in  Lee  et  al.  (1997),  which  is  attached  to  this  report  as  Appendix  1.  We  also  used  the 


same  control  algorithm  in  the  application  of  the  Lorentz  force  to  a  boundary  layer  control,  resulting 
in  a  substantial  improvement  over  an  open-loop  control  (Berger  et  al.,  2000). 

2.2  Suboptimal  Control 

In  another  approach  exploring  new  control  strategies,  we  developed  a  suboptimal  control  approach 
for  a  turbulent  channel  flow.  Instead  of  searching  for  a  global  optimal  condition,  for  which  time 
iterations  are  required  (and  hence  no  practical  importance),  the  suboptimal  approach  seek  an 
optimal  control  within  a  short  time  horizon,  thus  avoiding  time  iterations.  Two  simple  feedback 
laws  for  drag  reduction  were  obtained.  These  new  feedback  control  laws  required  pressure  or  shear- 
stress  information  at  the  wall.  More  practical  control  laws  requiring  only  the  local  distribution  of 
the  wall  pressure  or  one  component  of  the  wall  shear  stress  were  also  derived  and  were  shown  to 
work  equally  well.  Further  details  of  his  work  can  be  found  in  Lee  et  al.  (1998),  which  is  attached 
to  this  report  as  Appendix  2.  It  is  worth  mentioning  here  that  the  final  control  algorithm  obtained 
from  this  suboptimal  control  approach  was  very  similar  to  that  obtained  from  the  neural  network 
mentioned  above.  Both  control  schemes  reduced  the  viscous  drag  by  mitigating  the  strength  of 
near-wall  streamwise  vortices,  thus  confirming  our  earlier  premise  that  controlling  near-wall 
stremwise  vortices  is  a  key  to  viscous  drag  reduction  in  turbulent  boundary  layers. 


2.3  Systems  Theory  Approach 

A  third  approach  considered  was  a  system  theory  approach  to  flow  control.  The  main  objective  of 
this  approach  is  to  design  a  robust  controller  entirely  based  on  the  information  of  the  system  to  be 
controlled  rather  than  control  designer’s  physical  intuitions  or  prior  knowledge  of  the  system.  This 
approach  resulted  in  a  feedback  stabilization  to  delay  the  transition  in  a  laminar  channel  flow  to  a 
reduction  in  viscous  drag  in  a  fully  developed  turbulent  flow.  This  work  demonstrated  that 
although  the  system  theory  approach  is  based  on  a  linear  theory,  it  can  be  applied  to  fully  nonlinear 
flows.  Further  details  of  this  work  can  be  found  in  Joshi  et  al.  (1997),  which  is  attached  to  this 
report  as  Appendix  3.  In  Lee  et  al.  (2000),  a  reduced-order  linear  model  was  successfully  applied 
to  a  turbulent  channel  flow,  resulting  in  similar  viscous  drag  reduction  to  that  obtained  by  using  a 
neural  network  as  an  adaptive  controller  mentioned  earlier.  Motivated  by  the  success  of  the  linear 
system  theory  approach,  we  also  investigated  a  role  of  linear  processes  in  fully  nonlinear  flows.  It 
was  found  that  a  linear  process  plays  an  essential  role  in  maintaining  near-wall  turbulence  to  the 
extent  that  without  a  particular  linear  process,  turbulence  could  not  be  maintained  in  a  turbulent 
channel  flow.  Further  descriptions  of  this  work  can  be  found  in  Kim  and  Lim  (2000),  which  is 
attached  to  this  report  as  Appendix  4.  We  mention  in  passing  that  the  present  work  has  generated 
much  interest  in  the  application  of  control  theories  to  flow  problems  in  the  control  and  fluid 
mechanics  community  as  evidenced  by  several  recent  publications  on  the  subject  matter. 


2.4  Microflap  and  Immersed-Boundary  Method 

All  of  the  control  algorithms  mentioned  above  used  surface  blowing  and  suction  as  control  input. 
Although  such  actuation  is  simple  to  implement  in  numerical  experiments,  precise  control  of 
blowing  and  suction  in  practice  in  the  same  manner  as  implemented  in  the  numerical  experiments  is 


difficult  to  achieve.  In  order  to  develop  a  control  strategy  for  currently  available  actuators  such  as 
micro-flap  actuators,  development  of  a  new  numerical  approach,  known  as  an  immersed-boundary 
approach,  is  currently  underway.  A  brief  description  of  the  methdd  is  given  below. 

It  was  our  desire  to  simulate  flow  over  a  complex  boundary  while  retaining  the  simplicity  and 
efficiency  of  computation  in  a  Cartesian  system.  This  was  done  through  use  of  an  immersed 
boundary  method.  The  equations  of  fluid  motion  are  calculated  on  the  regular  geometry  of  a 
periodic  channel.  Insertion  of  a  microflap  into  the  flow  domain  creates  an  internal  boundary,  which 
must  be  simulated  properly.  We  used  an  idea  of  Peskin  and  McQueen  (1989)  of  using  a  body  force 
to  represent  the  dynamical  effect  on  the  fluid  of  a  moving  wall  immersed  in  the  fluid.  In  their 
method,  and  our  application,  the  position  of  the  internal  boundary  within  the  flow  domain  does  not 
generally  coincide  with  the  computational  mesh.  Flow  properties  are  interpolated  from  the 
computational  mesh  to  the  locations  of  boundary  nodes  in  order  to  calculate  the  force  of  the 
boundary  on  the  fluid.  Then  that  force  is  distributed  onto  nodes  of  the  computational  mesh  near  the 
internal  boundary. 

The  force  in  Peskin  and  McQueen's  work  was  based  on  a  balance  of  tension  in  an  elastic  wall 
against  the  hydrodynamic  forces.  We  instead  used  an  idea  developed  by  Mohd-Yusof  (1997)  of 
creating  a  fictitious  force  constructed  so  as  to  drive  the  velocity  at  desired  locations  to  desired 
values.  In  Mohd-Yusof  s  application,  the  desired  locations  are  a  riblet  surface  and  the  desired 
values  are  zero  to  fulfill  the  no-slip  condition.  In  our  application,  the  desired  locations  are  the  flap 
positions  and  the  desired  values  are  the  flap  velocities,  also  to  fulfill  the  no-slip  condition.  A 
difference  between  Mohd-Yusofs  application  and  ours  is  the  moving  boundary  created  by  the  flap. 
Another  difference  is  that  the  flap  is  an  internal  boundary  with  flow  on  both  sides.  Mohd-Yusofs 
riblet  surface  is  an  external  boundary  with  only  one  side  within  the  flow  domain. 

The  Navier-Stokes  equation  with  a  forcing  term  is  temporally  discretized  to  yield, 

i,n+1 

- f  u  •  Vu  =  —Vo  +vV2u  +  f.  (1) 

At 

In  a  neighborhood  of  the  flow  boundary,  we  already  know  what  we  want  the  velocity  to  be  at  the 
end  of  the  time  step. 


un+1  =  w,  (2) 

so  we  solve  for  the  body  force  that  would  bring  that  velocity  about 

f  =  — — —  +  u  •  Vu  +  S/p  -vV2u.  (3) 

At 

The  body  force  is  applied  at  two  sets  of  points,  those  just  below  the  immersed  boundary  and  those 
just  above.  Where  the  boundary  coincides  with  the  grid  the  body  force  is  applied  at  the  boundary 
and  at  the  points  just  below  it.  Application  of  the  force  at  additional  points,  particularly  below  the 
immersed  boundary,  was  tested  but  found  to  give  no  benefit.  Our  method  of  applying  the  body 


force  gives  flexibility  in  choosing  the  immersed  boundary  not  found  in  some  other  methods,  since 
there  is  no  need  to  line  up  the  boundary  with  a  grid. 

This  simulation  method  has  been  useful  for  examining  the  interaction  of  the  MEMS  actuator  with 
the  fluid.  Velocity  data  in  the  wind  tunnel  experiments  was  gathered  at  discrete  stations 
downstream  from  the  microflap.  Data  resolution  is  high  in  the  spanwise  and  wall-normal 
directions,  but  not  in  the  streamwise  direction.  The  numerical  simulations  make  data  at  any  point 
in  the  flow  available  for  interrogation.  The  ability  to  fill  in  the  gaps  of  the  wind  tunnel  data  has 
helped  explain  the  operation  of  the  actuator.  Wind  tunnel  data  showed  the  flap  causing  fluid  to 
move  away  from  the  wall  at  a  velocity  four  times  faster  than  the  flap  itself.  There  was  some 
question  as  to  whether  this  could  be  correct.  The  numerical  simulation  produced  similar  results  and 
provided  an  explanation.  Movement  of  the  flap  into  higher  speed  fluid,  combined  with  sucking 
fluid  under  the  opening  flap,  created  a  shear  layer,  which  would  roll  up  into  a  vortex  behind  the  flap 
and  eject  fluid  away  from  the  wall.  This  ejected  fluid  moved  away  from  the  wall  much  faster  than 
deflection  alone  could  produce.  This  observation  is  useful  for  a  couple  of  reasons.  First,  it  shows 
the  mechanisms  whereby  the  microflap  can  be  most  effective  in  influencing  the  flow.  Second,  it 
shows  that  a  flap's  actuation  will  be  felt  most  strongly  not  at  the  flap  itself,  but  instead  a  short 
distance  downstream. 

Laminar  flow  with  a  streamwise  vortex  passing  over  the  flap  was  generated  in  the  wind  tunnel 
experiments.  This  was  meant  to  be  a  simplified  analog  for  the  coherent  structures  found  near  the 
wall  in  a  boundary  layer.  The  vortex  pulls  high-speed  fluid  close  to  the  wall  and  generates  a  high 
shear  region  at  the  wall.  This  is  similar  to  the  sweep  events  in  turbulent  boundary  layers.  The 
microflap  actuator  counteracts  the  action  of  the  streamwise  vortex  and  reduces  shear  at  the  wall  in  a 
region  downstream  from  the  flap.  Our  numerical  simulations  give  results  similar  to  those  from  the 
wind  tunnel,  thus  adding  further  confidence  in  the  immersed-boundary  method  in  handling 
microflaps. 

Qualitative  comparisons  between  CFD  and  experimental  results  have  been  strong  and  encouraging. 
Quantitative  comparisons  are  nearing  completion,  and  when  successfully  completed  will  provide  a 
strong  validation  of  using  our  body  force  formulation  as  a  model  for  the  microflaps.  When  that 
work  is  completed,  simulation  efforts  will  turn  to  the  interaction  of  arrays  of  microflaps  with 
turbulent  channel  flows.  This  situation  has  not  yet  been  investigated  with  physical  experiments. 
Simpler  flap  models  will  also  be  compared  and  validated  against  the  body  force  formulation. 

Development  of  the  immersed-boundary  method  for  simulation  of  a  turbulent  flow  with  microflaps 
has  not  been  completed  yet.  However,  we  believe  that  the  immersed-boundary  approach  is  an 
efficient  and  cost-effective  approach  to  simulate  turbulent  flows  involving  complex  geometry.  We 
anticipate  that  the  immersed-boundary  method  will  emerge  as  a  powerful  tool  for  simulating  many 
turbulent  flows. 
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Application  of  neural  networks  to  turbulence  control  for  drag  reduction 
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A  new  adaptive  controller  based  on  a  neural  network  was  constructed  and  applied  to  turbulent 
channel  flow  for  drag  reduction.  A  simple  control  network,  which  employs  blowing  and  suction  at 
the  wall  based  only  on  the  wall-shear  stresses  in  the  spanwise  direction,  was  shown  to  reduce  the 
skin  friction  by  as  much  as  20%  in  direct  numerical  simulations  of  a  low-Reynolds  number  turbulent 
channel  flow.  Also,  a  stable  pattern  was  observed  in  the  distribution  of  weights  associated  with  the 
neural  network.  This  allowed  us  to  derive  a  simple  control  scheme  that  produced  the  same  amount 
of  drag  reduction.  This  simple  control  scheme  generates  optimum  wall  blowing  and  suction 
proportional  to  a  local  sum  of  the  wall-shear  stress  in  the  spanwise  direction.  The  distribution  of 
corresponding  weights  is  simple  and  localized,  thus  making  real  implementation  relatively  easy. 
Turbulence  characteristics  and  relevant  practical  issues  are  also  discussed.  ©  1997  American 
Institute  of  Physics.  [SI  070-663 1(97)02706-2] 


I.  INTRODUCTION 

The  ability  to  control  turbulent  flows  is  of  significant 
economic  interest.  Successful  control  of  turbulent  boundary 
layers  by  reducing  drag,  for  example,  can  result  in  a  substan¬ 
tial  reduction  in  operational  cost  for  commercial  aircraft  and 
marine  vehicles.  Recent  studies1,2  show  that  near- wall 
stream  wise  vortices  are  responsible  for  high  skin-friction 
drag  in  turbulent  boundary  layers.  Some  attempts  have  been 
made  to  reduce  the  skin-friction  drag  by  controlling  the  in¬ 
teractions  between  these  vortices  and  the  wall.  Choi  etal.,2 
for  example,  used  blowing  and  suction  at  the  wall  equal  and 
opposite  to  the  wall-normal  component  of  velocity  at 
y  + =  10.  They  showed  that  this  control  effectively  mitigated 
the  streamwise  vortices,  g;\  mg  approximately  25%  drag  re¬ 
duction  in  a  turbulent  channel  flow.  Although  the  method 
employed  in  their  work  is  impractical,  since  the  information 
aty+=  10  is  usually  not  available,  it  demonstrates  a  control 
scheme  which  reduces  skin-friction  drag  by  manipulation  of 
the  near-wall  streamwise  vortices.  Another  way  to  control 
the  streamwise  vortices  is  to  impose  spanwise  oscillation  by 
a  moving  wall  or  externally  imposed  body  force.3,4  These 
methods,  however,  require  a  large  amount  of  energy  input. 

A  systematic  approach  using  suboptimal  control  theory 
has  also  been  tried  in  the  past.  This  approach,  which  attempts 
to  minimize  a  cost  function,  was  successfully  applied  to  the 
stochastic  Burgers  equation.5  Moin  and  Bewley6  applied  a 
similar  approach  to  a  turbulent  channel  flow  to  achieve  up  to 
50%  drag  reduction.  The  control,  however,  requires  informa¬ 
tion  from  the  entire  velocity  field  inside  the  flow  domain  and 
excessive  computational  time,  making  it  impractical  to 
implement  in  real  situations.  For  practical  implementation,  a 
control  scheme  should  utilize  only  quantities  that  are  easily 
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measurable  at  the  wall,  and  should  be  fast  enough  to  be 
applied  in  real  time. 

The  objective  of  the  present  work  is  to  seek  wall  actua¬ 
tions,  in  the  form  of  blowing  and  suction  at  the  wall,  depen¬ 
dent  on  the  wall-shear  stress  to  achieve  a  substantial  skin- 
friction  reduction.  This  requires  knowledge  of  how  the  wall- 
shear  stresses  respond  to  wall  actuations,  i.e.,  the  correlation 
between  the  wall-shear  stresses  and  the  wall  actuations.  Be¬ 
cause  of  the  complexity  of  the  solutions  to  the  Navier- 
Stokes  equations,  .however,  it  is  not  possible  to  find  such  a 
correlation  in  closed  form  or  to  approximate  it  in  simple 
form.  Instead,  we  use  a  neural  network  to  approximate  the 
correlation  which  then  predicts  the  optimal  wall  actuations  to 
achieve  the  minimum  value  of  the  skin-friction  drag.  Neural 
networks  have  been  used  to  obtain  complicated,  nonlinear 
correlations  without  a  priori  knowledge  of  the  system  that  is 
to  be  controlled.  Jacobson  and  Reynolds,7  for  instance,  used 
a  neural  network  to  obtain  about  7%  drag  reduction  in  their 
simulation  of  an  artificial  flow.  In  this  paper,  we  describe 
how  we  constructed  and  trained  a  neural  network  off-line, 
and  then  implemented  an  on-line  control  scheme  (blowing 
and  suction  at  the  wall)  for  drag  reduction  based  on  that 
neural  network.  Applying  this  control  scheme  to  direct  nu¬ 
merical  simulations  of  turbulent  channel  flow  at  low  Rey¬ 
nolds  number,  we  observed  about  20%  drag  reduction.  We 
then  describe  how  examination  of  the  weight  distribution 
from  the  on-line  neural  network  led  to  a  very  simple  control 
scheme  that  worked  equally  well  while  being  computation¬ 
ally  more  efficient. 

In  Sec.  II,  a  brief  description  of  the  architecture  of  the 
neural  network  used  in  the  present  work  is  given.  In  Sec.  Ill, 
results  obtained  from  control  using  an  off-line  trained  net¬ 
work  are  presented,  while  results  obtained  from  control  using 
a  network  with  continuous  on-line  training  are  given  in  Sec. 
IV.  In  Sec.  V,  a  simple  control  law  based  on  the  weight 
distribution  from  the  successful  neural  network  control  is 
presented.  A  few  turbulence  statistics  are  given  in  Sec.  VI, 
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Surface  Shear  (dw/dy,) 
Stresses 

FIG.  1 .  Neural  network  architecture. 


followed  in  Sec.  VII  by  a  discussion  of  practical  implemen¬ 
tation  and  the  conclusions. 

In  this  paper  we  use  (x,y,z)  for  the  streamwise,  wall- 
normal,  and  span  wise  coordinates,  respectively,  and 
(w,u,w)  for  the  corresponding  velocity  components. 


II.  NEURAL  NETWORK 


In  this  section  we  describe  the  construction  of  a  neural 
network  to  learn  the  correlation  between  wall-shear  stresses 
and  wall  actuations  from  a  given  data  set.  Although  a  neural 
network  generally  requires  no  prior  knowledge  of  the  system 
(or  “plant”),  knowledge  about  the  near-wall  turbulence 
structures  provides  a  guideline  for  the  design  of  the  network 
architecture.  Initially  duldy  and  dw/dy  at  the  wall  at  several 
instances  of  time  were  used  as  input  data  fields  and  the  ac¬ 
tuation  at  the  wall  was  used  for  the  output  data  of  the  net¬ 
work.  Experimentally  we  found  that  only  dwldy  at  the  wall 
from  the  current  instance  of  time  was  necessary  for  sufficient 
network  performance.  Because  we  wanted  the  output  to  be 
based  only  on  a  local  input  area,  we  designed  our  network 
using  shared  weights.  The  network  had  a  single  set  of 
weights  (a  template)  that  is  convolved  over  the  entire  input 
space  to  generate  output  values;  that  is,  we  used  the  same  set 
of  weights  for  each  data  point  and  the  training  involves  iter¬ 
ating  over  all  data  points.  The  template  extracts  spatially 
invariant  correlations  between  input  and  output  data.  The 
size  of  the  template  was  initially  chosen  to  include  informa¬ 
tion  about  a  single  streak  and  streamwise  vortex,  and  then 
was  varied  to  find  an  optimal  size. 

We  used  a  standard  two-layer  feedforward  network  with 
hyperbolic  tangent  hidden  units  and  a  linear  output  unit  (see 
Fig.  1).  The  functional  form  of  our  final  neural  network  is 


{N—  1  )/2 


dw 


vik=Wa  tanh|  2  W,— 

/=-(  jv—  i  )/2  <*y 


-  w 

rr  c  , 


j,k+i 


l^k^Nz,  (1) 

where  the  W' s  denote  weights,  N  is  the  total  number  of  input 
weights,  and  the  subscripts  j  and  k  denote  the  numerical  grid 
point  at  the  wall  in  the  streamwise  and  spanwise  directions 
respectively.  Nx  and  N:  are  the  number  of  computational 


domain  grid  points  in  each  direction.  The  summation  is  done 
over  the  spanwise  direction.  Seven  neighboring  points 
(AT=  7),  including  the  point  of  interest,  in  the  spanwise  di¬ 
rection  (corresponding  to  approximately  90  wall  units  with 
our  numerical  resolution)  were  found  to  provide  enough  in¬ 
formation  to  adequately  train  and  control  the  near-wall  struc¬ 
tures  responsible  for  the  high-skin  friction.  Note  that  the 
blowing  and  suction  are  applied  at  each  grid  location  accord¬ 
ing  to  (1)  as  a  numerical  approximation  of  distributed  blow¬ 
ing  and  suction  on  the  surface.  A  scaled  conjugate  gradient 
learning  algorithm8  was  used  to  produce  rapid  training.  For 
given  pairs  of  (VjCk\dw/dy\jk),  network  was  trained  to  mini¬ 
mize  the  sum  of  a  weighted-squared  error  given  by 

error=  ^  2  ?)2,  (2) 

where  udcs  is  the  desired  output  value  and  unct  is  the  network 
output  value  given  by  Eq.  (1).  The  weights  were  initialized 
with  a  set  of  random  numbers.  Note  that  the  error  defined  in 
(2)  exponentially  emphasizes  (proportional  to  \)  large  actua¬ 
tions.  This  error  scaling  was  chosen  based  on  Choi  et  aV  s~ 
observation  that  large  actuations  are  more  important  for  drag 
reduction.  Usually  within  100  training  epochs,  the  error 
reached  its  asymptotic  limit. 


III.  OFF-LINE  TRAINING  AND  CONTROL 

As  an  initial  experiment  we  investigated  whether  we 
could  train  a  neural  network  to  predict  the  velocity  at 
y+=  10  from  only  the  wall-shear  stresses.  The  rationale  be¬ 
hind  this  experiment  was  to  duplicate  an  existing  control  law 
using  only  measurements  available  at  the  wall  such  that  the 
output  from  the  network  could  be  used  as  input  to  the  actua¬ 
tor.  The  network  should  yield  a  similar  amount  of  drag  re¬ 
duction  to  that  obtained  by  Choi  et  al 2  through  prediction  of 
the  velocity  at  y+  =  10.  The  training  data  consisted  of  100 
time  steps  of  output  obtained  from  a  numerical  simulation  of 
channel  flow  under  the  control  employed  by  Choi  et  al .,2  i.e., 
using  the  wall-normal  velocity  at  y  +  =  10.  The  flow  regime  is 
turbulent  channel  flow  with  Rer=  100,  where  ReT  is  the  Rey¬ 
nolds  number  based  on  the  wall-shear  velocity,  uT ,  and  the 
channel  half-width,  S.  All  numerical  simulations  presented  in 
this  paper  were  obtained  using  a  modified  version  of  Kim 
et  al.  ’s9  spectral  code  with  the  computational  domain 
(4  7r, 2,4  7t/3)<5,  and  a  grid  resolution  of  (32,  65,  32)  in  the 
(x9y,z)  directions,  respectively.  Each  time  step  contained  a 
32X32  array  of  input  values  (dwldy\w)  and  corresponding 
actuations,  —  v  aty  +  =10. 

We  trained  several  networks  with  one  hidden  unit  and 
different  sized  input  templates:  7X1,  7X3,  7X5,  9X1, 
9x3,  1 1 X  1 ,  and  11X3  (the  number  of  input  points  in  the 
spanwise  direction  by  the  number  of  input  points  in  the 
streamwise  direction).  After  training  was  completed,  the  dis¬ 
tribution  of  input  weights  was  examined  to  see  whether  there 
was  a  discernible  pattern  in  the  input  template.  The  weight 
distributions  in  the  spanwise  direction  at  the  same  stream- 
wise  location  for  different  input  templates  are  shown  in  Fig. 
2.  The  same  pattern  for  all  7  input  template  sizes  was  ob¬ 
served,  and  is  similar  to  a  finite  difference  approximation  of 
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FIG.  2.  Weight  distribution  from  off-line  trainings  for  various  sizes  of  the 
input  template:  7  X  1,7  X  3,7  X  5,9  X  1,9  X  3, 11  X  1,11  X  3.Weightsare 
normalized  by  Wx. 


spanwise  differentiation.  Increasing  the  number  of  hidden 
units  only  marginally  improved  performance,  whereas  in¬ 
creasing  the  template  size  significantly  reduced  the  final 
training  error. 

We  then  applied  a  control  scheme  to  a  regular  channel 
flow  by  fixing  the  final  input  weights  obtained  from  the  off¬ 
line  training.  This  was  perhaps  a  somewhat  naive  approach 
since  the  weights  were  obtained  from  fully  controlled  flow 
data  and  the  shear  stresses  used  for  the  training  were  already 
altered  by  the  actuations.  Nevertheless,  two  cases  were 
tested:  a  control  scheme  based  on  7  weights  in  the  spanwise 
direction,  and  another  based  on  the  same  7  weights  plus  3 
more  in  the  immediate  downstream  location.  These  10 
weights  were  chosen  because  among  all  weights  obtained 
from  the  off-line  training  they  had  non-negligible  values.  In 
Fig.  3,  the  mean  shear  stress  variations  at  the  wall  (i.e.,  drag) 
obtained  with  these  two  controls  are  plotted  along  with  the 
no-control  case.  Nearly  18%  drag  reduction  was  achieved, 
with  slightly  better  performance  from  the  10-weight  net¬ 
work. 

These  results  demonstrate  that  a  correlation  exists  be¬ 
tween  the  shear  stresses  at  the  wall  and  the  desired  actua- 


FIG.  3.  Mean  wall -shear  stress  histories  for  various  control  laws  compared 
to  the  no-control  case:  — ,  no  control;  — ,  control  with  7  fixed  weights 
obtained  from  off-line  control;  •••,  control  with  10  fixed  weights  from  the 
same  off-line  control.  All  runs  with  control  shown  in  this  paper  were  inte¬ 
grated  over  at  least  20  time  units  that  is  much  longer  than  the  typical  turn¬ 
over  time  of  the  streamwise  vortices,  1  time  unit. 


FIG.  4.  Schematic  representation  of  adaptive  inverse  model  control. 


tions,  and  that  control  based  on  this  correlation  produces  a 
significant  amount  of  drag  reduction.  This  fixed-weight  con¬ 
trol  scheme,  however,  was  deduced  from  fully  controlled 
flow  data;  thus  it  does  not  guarantee  the  same  performance 
for  other  flow  situations. 


IV.  ON-LINE  CONTROL 


In  the  previous  section  we  showed  how  successful  con¬ 
trol  based  on  off-line  training  can  be  obtained.  However, 
since  the  system  we  are  trying  to  control  is  time- varying  and 
nonlinear,  this  approach  is  not  likely  to  generalize  well.  Con¬ 
tinuous  on-line  training  allows  a  controller  to  adapt  to  the 
evolution  of  the  system.  In  this  section,  we  describe  an  adap¬ 
tive  controller  for  on-line  training  and  control. 

There  are  various  schemes  for  on-line  neural  network 
control.  The  most  direct  scheme  is  adaptive-inverse  model 
control.10  A  schematic  representation  of  this  approach  is 
shown  in  Fig.  4.  Here  the  plant  denotes  the  numerical  solver 
of  the  Navier-Stokes  equations.  This  configuration  employs 
a  neural  network  to  model  the  (possibly  time-varying)  in¬ 
verse  plant  mapping  from  wall-shear  stress  dwldy  |H,  to  the 
wall-normal  actuations,  and  then  uses  a  copy  of  the  model  as 
the  controller  with  thw  desired  shear  stresses  as  input.  One 
restriction  of  this  technique  is  that  it  usually  requires  an  ini¬ 
tial  model  training  phase  using  random  plant  inputs  and  cor¬ 
responding  plant  outputs.  This,  however,  caused  no  serious 
problems  since  usually  one  time  step  was  enough  for  model 
training  because  of  the  shared-weight  architecture  in  our  ap¬ 
plication  (each  time  is  1024  data  points  and  our  networks 
only  have  a  few  weights).  Once  the  model  represents  a  rea¬ 
sonably  close  approximation  to  the  actual  plant  inverse,  a 
copy  is  then  implemented  as  a  feedforward  controller. 

The  desired  inputs  to  the  controller  are  a  fractional  re¬ 
duction  in  the  shear  stress  from  the  previous  step,  i.e., 


i  dw\des  I  dw\ 

\  dy  /  \  dy] 

v  s  >  i  +  A/  '  7  '  t 


(3) 


where  0<  rj<  1 .  Based  on  the  networks  and  techniques  we 
investigated,  this  indirect  suppression  of  dwldy  |w,  instead  of 
duldy |H. ,  turns  out  to  be  more  efficient  in  achieving  drag 
reduction.  The  output  of  the  controller,  which  is  the  input  to 
the  plant,  is  the  predicted  actuation  necessary  to  produce  this 
shear-stress  reduction.  The  quantity  77  should  be  chosen  such 
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FIG.  5.  Distribution  of  the  weights:  0 ,  from  on-line  training;  — , 
A(  1  -cos(7T/))//.  Az  +  =  13. 


that  sets  of  the  desired  large  amplitude  outputs  are  among  the 
sets  that  are  well  represented  by  the  training  sets.  Good  per¬ 
formance  was  achieved  for  the  range  of  77=0.8—0.85. 

A  turbulent  channel  flow  at  Reynolds  number 
ReT=  100  was  used  to  test  the  neural  network.  We  allowed 
all  the  weights  in  the  network  to  adapt  and  examined  the 
input  template  pattern  after  each  time  step.  As  the  control 
began,  the  weight  distribution  immediately  assumed  a  fixed 
pattern  similar  to  the  off-line  pattern.  There  was  no  appre¬ 
ciable  change  in  the  relative  magnitudes  of  the  template 
weights  over  time,  indicating  that  the  pattern  is  preserved. 
The  absolute  magnitudes,  however,  did  vary  indicating  the 
need  for  gain  and  bias  adaptation  for  each  layer.  The  number 
of  hidden  layers,  the  number  of  the  hidden  units,  the  size  of 
the  input  template,  and  the  error  scale  of  the  training  error 
[X  of  Eq.  (2)]  were  1,  1,  7X1,  and  5,  respectively.  We 
varied  the  values  of  the  error  scale  and  found  5  to  be  opti¬ 
mum,  with  larger  values  causing  an  instability  in  training. 
We  also  varied  the  number  of  the  hidden  units,  but  the  net¬ 
work  with  a  single  hidden  unit  produced  the  best  result.  The 
converged  weight  distributions  for  20  consecutive  time  steps 
after  an  asymptotic  state  was  reached  are  shown  in  Fig.  5. 
The  pattern  is  only  slightly  different  from  the  one  obtained 
from  the  off-line  training  (see  Fig.  2).  It  should  be  noted  that 
this  pattern  emerged  immediately  after  the  on-line  control 
began. 

Time  histories  of  the  wall-shear  stress  for  3  different 
input  template  sizes  are  shown  in  Fig.  6.  All  other  network 
parameters  are  kept  the  same.  The  abscissa  represents  non- 
dimensional  time  normalized  by  uT  and  S.  One  time  unit 
roughly  corresponds  to  the  time  for  a  fluid  particle  traveling 
with  the  centerline  velocity  to  move  about  20  channel  half 
widths.  The  computations  were  carried  out  long  enough  to 
ensure  that  a  statistically  steady  state  was  reached.  As  the 
control  began,  the  drag  quickly  drops  to  about  80%  of  that 
observed  without  control,  for  the  two  cases  with  template 
sizes  of  7  X  1  and  9X1.  The  template  size  of  5  X  1 ,  however, 
did  not  produce  as  much  reduction,  indicating  that  at  least  7 
spanwise  points  (which  extended  about  90  wall  units,  with 
our  grid  resolution)  should  be  used  for  good  performance.  At 
the  initial  stage  of  our  study,  we  tested  a  control  using  both 
dxv/dy |H.  and  duldy\„  as  input,  with  a  7X5  input  template. 
It  produced  about  the  same  amount  of  reduction,  but  with 
excessive  training  time  due  to  a  significantly  larger  number 


FIG.  6.  Mean  wall-shear  stress  histories  for  on-line  controPj^th  different 
input  template  sizes:  — ,  no  control;  control  with  5X1;  — ,  control  with 
7X1;  — ,  control  with  9X  1. 


of  weights.  Furthermore,  it  did  not  produce  coherent  patterns 
in  the  weight  distributions. 

Since  the  7  weights  showed  a  fixed  pattern,  we  fixed  the 
input  template  weights  to  this  pattern  and  used  a  single  hid¬ 
den  unit  network,  giving  only  4  adaptable  parameters  (a  bias 
and  gain  for  each  layer).  This  simplified  network  had  the 
following  functional  form; 

vjt=Watznh(Wbg-Wc)-Wd  (4) 

with 


^  <?w 


j.k+i 


(5) 


where  W^s  are  the  fixed-weight  pattern  obtained  from  the 
previous  on-line  control.  On-line  control  using  this  network 
produced  a  similar  amount  of  drag  reduction.  The  weight 
variations  with  time  were  also  monitored.  The  bias  weights 
( Wc  and  W d)  were  negligibly  small.  Those  controlling  the 
gain  (Wa  and  Wb)y  which  had  finite  values,  changed  in. time 
significantly,  although  the  product  of  the  two  gaW  -Weights 
remained  almost  constant.  This  suggests  that  effective  con¬ 
trol  can  be  achieved  by  simply  using  g  with  an  adjustable 
amplitude.  This  will  be  discussed  in  the  following  section. 


V.  A  SIMPLE  CONTROL  SCHEME 


Since  the  control  based  on  the  network  given  by  Eq.  (4) 
produced  a  substantial  reduction  in  drag,  this  section  devel¬ 
ops  a  control  scheme  based  only  on  the  weighted  sum  of  the 
wall-shear  stress.  The  distribution  of  the  weights  can  be  ap¬ 
proximated  by  (see  Fig.  5): 


1  -COS(777) 

W*=A - : - 

;  J 


(6) 


where  j~0  corresponds  to  the  point  where  the  control  is 
applied.  Since  only  the  relative  values  are  important,  the  con¬ 
stant  A  has  no  special  meaning.  A  computation  with  a  twice 
higher  resolution  was  run  to  confirm  that  Eq.  (6)  gives  the 
proper  form  of  the  weight  distribution  (see  Fig.  7).  It  turns 
out  that  Eq.  (6)  is  the  inverse  Fourier  transform  of  ikzl\kz\, 
where  kz  is  the  wave  number  in  the  spanwise  direction,  for 
the  finite  maximum  wave  number,  i.e., 
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FIG.  7.  Distribution  of  the  weights:  0,  from  on-line  training;— 
A  ( 1  -  cos(7J j))fj.  Az  +  -  6 . 


140 

130 

120 

$L  no 
100 

90 

80 

70 

202530  35  40455055 

t 

FIG.  9.  Mean  wall-shear  stress  histories  for  various  control  schemes  com¬ 
pared  to  the  no-control  case:  — ,  no  control;  — ,  on-line  control  with  neural 
network  with  7X  1  template;  control  with  7  fixed  weights. 


r*m  i kz 

J-JO 


exp(  —  ikzz)dkz  =  2 


1-co  s(kmz) 
z 


(7) 


where  km-  tt l Az  is  the  maximum  wave  number  and  A z  is 
the  numerical  grid  spacing.  Replacing  z  with  /Az  in  the 
right-hand-side  of  Eq.  (7)  leads  to  Eq.  (6).  From  this  result 
and  the  convolution  theorem,  one  can  suggest  the  following 
simple  control  law  for  wall  transpiration: 


ikz  dw 

Vw=C\ t\~*y 


1  d 

/  dw\ 

|*.|  dz 

H-  1  *  1 

\*yl 

(8) 


where  the  “hat”  denotes  a  Fourier  transformed  quantity  and 
C  is  a  positive  scale  factor  determining  the  amplitude  of  the 
actuation.  Equation  (8),  which  should  produce  similar  drag 
reduction  as  the  result  with  7  fixed  weights  [Eq.  (4)],  implies 
that  the  optimum  blowing  and  suction  at  the  wall  is  propor¬ 
tional  to  d(dw!dy)/dz  with  the  high  wave  number  compo¬ 
nent  suppressed  by  1/|A:-|.  Note  that  the  wall-normal  actua¬ 
tion  proportional  to  d{dw!  dy)!  dz  counteracts  the  up-and- 
down  motion  induced  by  a  streamwise  vortex  (see  Fig.  8), 
consistent  with  that  used  by  Choi  et  al 2  Smith  et  alu  dis¬ 
cussed  the  stability  of  the  near-wall  turbulence  structure  to  a 
local  adverse  pressure  gradient  imposed  near  the  surface. 


FIG.  8.  A  schematic  of  a  flow  field  induced  by  a  streamwise  vortex  and 
corresponding  d{dw!  dy\,J)l  dz. 


The  suction  at  the  left  side  of  the  vortex  in  Fig.  8  can  relieve 
such  adverse  pressure  gradient,  thus  preventing  the  surface- 
layer  eruption.  Choi  et  al 2  reported  that  actuation  based  only 
on  the  local  value  of  8{dw! dy)!  dz,  which  is  also  a  dominant 
term  of  the  Taylor  expansion  of  v  at  some  distance  away 
from  the  wall,  did  not  produce  such  a  substantial  drag  reduc¬ 
tion.  This  indicates  that  suppression  of  high  wave  number 
components  by  V\kz\9  or  equivalently  using  locally  weighted 
value  of  dw/dy,  plays  a  key  role  in  reducing  drag.  Note  that 
the  weights  at  even  numbered  grid  points  away  from  the 
center  point  vanish.  This  is  beneficial  for  physical  implemen¬ 
tation,  since  sensors  and  actuators  cannot  be  placed  at  the 
same  location.12  Because  the  Fourier  integral  is  computed  for 
a  finite  value  of  km ,  the  values  of  the  weight  at  noninteger 
j  in  Eq.  (6)  have  no  meaning.  Equation  (8)  is  equivalent  to 


<^,/2  dw 
Vjk=C  2  Wt— 

i—  —  (N—  1  )/2  oy 


j.k+i 


(9) 


where  Wt  is  given  by  Eq.  (6).  The  magnitude  of  the  weights 
decays  with  increasing  distance  from  the  center,  which  al¬ 
lows  for  good  approximation  using  only  a  small  number  of 
weights,  i.e.,  successful  control  requires  only  local  values  of 
the  shear  stress.  A  natural  concern  is  how  different  grid  spac- 
ings  affect  the  control.  Since  the  above  control  law  [Eq.  (9)] 
is  simply  another  expression  of  Eq.  (8),  which  is  a  good 
approximation  as  long  as  km  is  large  enough,  the  control 
should  be  relatively  independent  of  resolution.  This  was  con¬ 
firmed  by  a  computation  with  a  higher  resolution. 

Control  based  on  Eq.  (9)  with  7  points  produced  the 
same  amount  of  drag  reduction  (20%)  as  the  neural  network 
control  (see  Fig.  9).  The  constant  C  is  chosen  so  that  the 
root-mean-squared  (rms)  value  of  the  actuation  is  kept  at 
0.15kt.  Blowing  and  suction  of  this  magnitude  at  the  wall 
suppress  the  near-wall  streamwise  vortices  by  counteracting 
the  up-and-down  motions  associated  with  these  vortices  (see 
Fig.  8). 

The  fact  that  a  control  scheme  that  is  essentially  linear 
[Eq.  (9)]  produced  the  same  amount  of  drag  reduction  as 
on-line  control  with  a  nonlinear  network  led  to  the  following 
question:  Can  the  same  weight  pattern  be  captured  when 
on-line  control  with  a  “linear”  network  is  used?  To  answer 
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FIG.  10.  Contours  of  streamwise  vorticity  in  a  cross-flow  plane:  (a)  no 
control;  (b)  control  using  7  fixed  weights.  The  contour  level  increment  is  the 
same  for  both  figures.  Negative  contours  are  chain-dotted. 


this  question,  we  tested  an  on-line  controller  utilizing  a  linear 
network.  The  hyperbolic  tangent  function  in  Eq.  (1)  was  re¬ 
placed  with  a  linear  function  and  the  hidden  layer  was  re¬ 
moved.  All  other  parameters  were  kept  the  same  as  the  non¬ 
linear  network  case.  The  linear  network  seemed  to  produce 
almost  the  same  pattern  in  weights  as  the  one  produced  by 
our  nonlinear  network.  Also,  drag  was  reduced  by  the  same 
amount.  The  pattern  that  was  obtained  by  averaging  over 
time,  however,  was  not  well  preserved  in  time.  The  weight 
W3  of  Eq.  (1),  for  example,  frequently  showed  order  of  mag¬ 
nitude  excursions  resultinp  in  a  large  standard  deviation,  1.5, 
compared  to  the  value  of  the  weight,  0.3  (weights  were  nor¬ 
malized  by  the  maximum  amplitude  fVj).  The  nonlinear  net¬ 
work,  however,  produced  more  constant  weights;  the  stan¬ 
dard  deviation  of  W3  was  only  0.15.  This  result  indicates  that 
the  nonlinearity  should  add  robustness  by  bounding  the 
weights.  This  can  simplify  hardware  implementation  by  lim¬ 
iting  the  signals  and  subsequently  the  necessary  dynamic 
range  of  the  processing  circuitry. 


VI.  TURBULENCE  STATISTICS 

The  computed  flow  fields  for  a  no-control  case  and  a 
successful  control  case,  based  on  the  7-point  weighted  sum 
of  Sw/dy  |H.  [Eq.  (9)],  were  examined  to  investigate  the 
mechanism  by  which  the  drag  reduction  is  achieved.  The' 
most  salient  feature  of  the  controlled  case  was  that  the 
strength  of  the  near-wall  streamwise  vortices  was  drastically 
reduced.  In  Fig.  10,  contours  of  streamwise  vorticity  in  a 
cross  plane  are  shown.  This  result  further  substantiates  the 
notion  tint  a  successful  suppression  of  the  near-wall  stream- 
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FIG.  1 1.  Probability-density  function  of  the  wall-shear  stress:  — ,  no  con¬ 
trol;  — ,  control  with  7  fixed  weights.  The  area  under  each  curve  is  normal¬ 
ized  to  one. 

wise  vortices  leads  to  a  significant  reduction  in  drag.  Note 
that  for  the  controlled  case  the  wall  actuations  were  applied 
at  both  walls. 

The  probability-density  function  of  the  wall-shear  stress 
in  the  streamwise  direction  is  shown  in  Fig.  11.  It  is  evident 
that  the  control  is  very  effective  in  suppressing  large  fluctua¬ 
tions,  thus  reducing  the  mean-skin  friction.  Furthermore,  the 
root-mean-square  (rms)  values  of  turbulent  fluctuations  in 
the  wall  region  are  also  reduced  as  shown  in  Fig.  12(a).  The 
same  trend  was  observed  by  Choi  et  al. 2  The  rms  vorticity 
fluctuations  are  also  significantly  reduced,  except  for  ojx  very 


(a) 


FIG.  12.  Root-mean-square  fluctuations  normalized  by  the  wall  variables: 
— ,  no  control;  control  with  7  fixed  weights,  (a)  Velocity  fluctuations;  (b) 

vorticity  fluctuations. 
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FIG.  13.  Contours  of  the  wall  actuation:  (a)  control  using  dw!dy\w  with  7 
fixed  weights;  (b)  control  using  information  at  y+ - 10.  Negative  contours 
are  chain-dotted. 

close  to  the  wall.  The  increase  for  the  latter  is  caused  by 
additional  dvldz  due  to  the  wall  actuation.  The  rms  fluctua¬ 
tions  of  <x)z  at  the  wall,  which  is  mainly  due  to  duldy  at  the 
wall,  is  also  decreased.  This  indicates  that  the  control  scheme 
led  to  a  reduction  in  the  mean  shear  and  its  variance  at  the 
wall  by  suppressing  large  fluctuations  as  shown  in  Fig.  1 1 . 
The  reduction  in  the  rms  fluctuations  in  a)x  and  co}  indicates 
that  the  control  scheme  indeed  reduced  the  strength  of  the 
near-wall  streamwise  vortices  and  the  wall-layer  streaks. 

Figure  13  compares  the  distribution  of  wall  actor  don 
used  in  our  control  with  that  from  Choi  et  al.'s2  u-control 
using  the  information  at  j/+=10  for  the  same  wall-shear 
stress  distribution.  The  corresponding  wall-shear  stress  dis¬ 
tribution  is  also  shown  in  Fig.  13.  The  wall  actuations  indi¬ 
cate  a  strikingly  similar  distribution  to  each  other,  even 
though  the  wall  actuation  of  our  control  is  based  only  on  the 
wall-shear  stress  dw! dy  |  H. .  Basically  our  control  scheme  de¬ 
tects  edges  of  locally  high-shear  stress  regions  by  measuring 
the  spanwise  variation  of  dwldy,  and  applies  appropriate 
wall  actuation  as  shown  in  Figs.  13(a)  and  13(b).  Since  high- 
shear  stress  regions  are  usually  elongated  in  the  streamwise 
direction,  only  spanwise  variation  is  necessary  for  detecting 
the  edges.  Since  dwldy  is  a  direct  measurement  of  stream- 
wise  vortices  (see  Fig.  8),  it  provides  more  appropriate  infor¬ 
mation  than  duldy.  This  is  consistent  with  our  finding  that 
dwldy  at  several  points  in  spanwise  direction  are  enough  for 
good  performance  of  control. 

VIL  DISCUSSION  AND  CONCLUSION 

We  have  presented  a  successful  application  of  a  neural 
network  to  turbulence  control  for  drag  reduction.  First  we 


were  able  to  construct  and  train  a  neural  network  off-line  to 
find  a  correlation  between  the  wall-shear  stress  and  the  de¬ 
sired  wall  actuations  coming  from  the  information  at 
y+=  10.  Based  on  the  optimal  network  structure  from  the 
off-line  training,  we  successfully  implemented  an  on-line  in¬ 
verse  model  controller  in  numerical  experiments  of  a  turbu¬ 
lent  channel  flow,  resulting  in  about  20%  drag  reduction. 
Finally,  we  were  able  to  deduce  a  simple  control  scheme 
[Eq.  (9)]  for  drag  reduction  based  on  observation  of  .the 
weight  distribution  from  a  successful  control  case.  This  con¬ 
trol  scheme  uses  a  minimal  amount  of  wall -shear  stress  in¬ 
formation  and  requires  only  simple  operations,  thus  render¬ 
ing  a  scheme  whose  actual  hardware  implementation  would 
be  relatively  easy. 

Beginning  with  a  nonlinear  network,  we  found  a  simple 
linear  control  scheme  from  the  pattern  of  the  weights,  sug¬ 
gesting  the  use  of  a  linear  network.  Although  a  linear  net¬ 
work  produced  a  similar  pattern,  the  pattern  was  not  well 
preserved  in  time.  The  weights  can  vary  significantly  in¬ 
creasing  the  difficulty  in  circuit  implementation  due  to  the 
larger  required  dynamic  range.  The  nonlinear  network,  how¬ 
ever,  produced  bounded  weights  simplifying  hardware 
implementation.  This  suggests  that  a  certain  amount  of  non¬ 
linearity  is  beneficial  to  capture  a  stable  coherent  pattern  for 
our  system. 

Although  the  present  control  scheme  is  a  significant  im¬ 
provement  over  earlier  approaches  that  require  velocity  in¬ 
formation  within  the  flow,  there  is  a  technical  issue  before 
the  present  control  scheme  can  be  implemented  in  real  prac¬ 
tice.  Precise  control  of  blowing  and  suction  distributed  over 
a  surface  is  difficult  to  implement.  Other  approaches,  such  as 
surface  movements  by  deformable  surface,  may  prove  to  be 
more  practical. 

There  are  also  several  fundamental  issues  that  must  be 
addressed.  First,  our  numerical  experiments  were  performed 
in  a  very  low  Reynolds  number  flow.  It  remains  to  be  seen 
whether  the  same  control  scheme  extends  to  higher  Reynolds 
numbers.  If  the  main  cause  of  high-skin  friction  in  turbulent 
boundary  layers  at  higher  Reynolds  numbers  is  also  due  to 
the  near-wall  streamwise  vortices,  the  same  scheme  should 
work  equally  well.  Detection  of  these  vortices  through  the 
wall-shear  stress  would  become  more  difficult,  however,  be¬ 
cause  the  scales  associated  with  these  vortices  decrease  as 
the  Reynolds  number  increases. 

Another  important  issue  is  the  influence  of  the  spatial 
resolution  of  sensors  and  actuators  on  performance.  We 
showed  that  a  simple  control  scheme  with  4  neighboring 
points  in  the  spanwise  direction  (3  among  7  weights  are  ze¬ 
ros),  which  span  90  wall  units,  performed  very  well.  This 
suggests  that  the  distance  between  sensors  should  be  about 
20  wall  units.  As  the  Reynolds  number  increases,  the  physi¬ 
cal  separation  between  sensors  must  decrease.  We  note  in 
passing  that  recent  advances  in  micromachined  sensors  and 
actuators  make  these  scales  feasible,13  and  the  present  work 
is  part  of  a  joint  research  project  aimed  at  integrating  micron¬ 
sized  sensors  and  actuators  with  analog  VLSI  control  cir¬ 
cuitry. 

A  third  issue  is  determining  an  appropriate  scale  factor 
C  in  Eq.  (9).  We  tested  a  range  of  C  and  empirically  found 
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that  values  yielding  an  actuation  rms  value  between  0.1  ur 
and  0.15wt  gave  the  best  performance.  Smaller  values  re¬ 
sulted  in  less  reduction,  while  larger  ones  caused  rapid  fluc¬ 
tuations  of  the  wall-shear  stress  over  time.  The  experiments 
were  only  performed  for  Rer=100  and  180,  from  which  the 
optimum  value  was  deduced,  but  we  expect  that  the  same 
amplitude  level  should  produce  similar  reductions  for  higher 
Reynolds  number  flows.  We  further  note  that  for  this  ampli¬ 
tude  the  required  power  input  per  unit  area  to  produce  the 
actuation,  p„vw+0.5pvl{~0.lpu3T )  was  significantly  less 
compared  to  the  power  saved  due  to  the  reduced  drag, 
bCfICf  rwUc(~~3.5pul)  for  20%  reduction  in  Cf .  Here 
pw,  p,  Cf ,  rH.  and  Uc  are  the  wall  pressure,  density,  friction 
coefficient,  averaged  wall-shear  stress  in  the  streamwise  di¬ 
rection,  and  centerline  velocity,  respectively.  Note  that  we 
did  not  account  for  device  loss  to  deliver  blowing  and  suc¬ 
tion  nor  the  power  consumed  by  sensors  in  this  estimate. 

One  final  practical  issue  worth  mentioning  is  the  time 
delay  between  sensing  and  actuation.  None  was  included  in 
any  of  our  numerical  experiments.  In  a  real  situation,  how¬ 
ever,  there  will  be  a  finite  time  delay  between  sensing  and 
actuation.  The  process  time  for  the  simple  control  scheme 
should  be  negligible  since  a  weighted  summation  is  an  in¬ 
volved  process,  thus  not  imposing  any  restriction  on  the  time 
delay.  Ideally,  sensors  should  register  the  response  of  the 
near-wall  structures  to  be  controlled  due  to  wall  actuations, 
which  will  take  a  certain  amount  of  time.  Blowing  and  suc¬ 
tion  at  the  wall,  however,  can  produce  an  immediate  spurious 
response  at  nearby  sensors  due  to  local  conservation.  One 
possible  remedy  to  this  problem  is  an  under-relaxation  of  the 
actuation  signal  with  past  signals,  accounting  for  the  re¬ 
sponse  of  a  longer  time  scale.  For  example,  an  underrelax¬ 
ation  scheme  using  the  following  formula 


w( 


Bw 

By 


+  (1-£K 


JM 


jk 


(10) 


could  take  into  account  all  past  signals  with  a  single  param¬ 
eter  £  (0<£<  1 ).  We  used  this  approach  successfully  in  our 
numerical  experiment  for  Rer=  180  and  found  that  there  is 
an  optimum  £  depending  on  temporal  resolution. 

The  most  significant  finding  of  the  present  study  is  that  a 
single  span  wise  strip  of  Bwldy  |H,  is  enough  to  achieve  sig¬ 
nificant  drag  reduction.  Additional  information  about 
dwldy\y.  in  the  streamwise  direction  or  Bui  By  |H.  in  either 
direction  reduced  the  efficiency  of  our  neural-network  based 
control.  We  also  tested  different-sized  input  templates  and 
found  that  as  the  template  size  increased  in  the  streamwise 
direction,  the  training  time  became  excessive  and  the  fixed 
pattern  in  the  weight  disappeared,  indicating  that  too  many 
unnecessary  weights  caused  training  problems. 
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Two  simple  feedback  control  laws  for  drag  reduction  are  derived  by  applying  a 
suboptimal  control  theory  to  a  turbulent  channel  flow.  These  new  feedback  control 
laws  require  pressure  or  shear-stress  information  only  at  the  wall,  and  when  applied 
to  a  turbulent  channel  flow  at  Rer  =  110,  they  result  in  16-22%  reduction  in  the 
skin-friction  drag.  More  practical  control  laws  requiring  only  the  local  distribution  of 
the  wall  pressure  or  one  component  of  the  wall  shear  stress  are  also  derived  and  are 
shown  to  work  equally  well. 


1.  Introduction 

Recent  studies  have  shown  that  near-wall  streamwise  vortices  are  responsible 
for  high  skin-friction  drag  in  turbulent  boundary  layers.  Many  attempts  aiming  at 
controlling  these  vortices  have  been  made  in  order  to  achieve  a  skin-friction  drag 
reduction  in  turbulent  boundary  layers.  Most  of  such  attempts,  however,  have  been 
ad  hoc,  largely  based  on  physical  intuition.  The  active  control  by  Choi,  Moin  &  Kim 
(1994),  for  example,  used  blowing  and  suction  at  the  wall  that  is  equal  and  opposite 
to  the  wall-normal  component  of  the  velocity  at  >,+  =  10,  resulting  in  as  much  as 
25%  reduction  in  their  numerical  simulation.  This  approach,  however,  is  impractical 
since  the  required  velocity  information  at  y+  =  10  is  not  normally  available.  For 
any  practical  implementation  a  control  scheme  should  be  based  solely  on  quantities 
measurable  at  the  wall. 

A  more  systematic  approach  based  on  an  optimal  control  theory  was  proposed 
by  Abergel  &  Temam  (1990).  Choi  et  al  (1993)  proposed  a  ‘suboptimal’  control 
procedure,  in  which  the  iterations  required  for  a  global  optimal  control  were  avoided 
by  seeking  an  optimal  condition  over  a  short  time  period.  The  suboptimal  control 
procedure  was  successfully  applied  to  control  of  the  Burgers  equation.  Bewley  & 
Moin  (1994)  were  the  first  to  apply  the  suboptimal  control  procedure  to  a  turbulent 
flow  and  reported  about  17%  drag  reduction.  The  procedure  developed  by  Bewley 
&  Moin  (1994)  still  requires  velocity  information  inside  the  flow  in  order  to  solve 
the  adjoint  problem,  from  which  a  feedback  control  input  was  derived.  In  spite  of 
this  obvious  drawback,  however,  the  fact  that  a  control  theory  applied  to  a  turbulent 
flow  resulted  in  a  substantial  drag  reduction  is  encouraging,  since  their  control 
procedure  was  derived  rigorously  from  a  control  theory,  in  which  a  pre-determined 

t  Present  address:  Department  of  Mechanical  Engineering,  University  of  Seoul,  Seoul,  Korea. 
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cost  functional  was  minimized.  Hill  (1993,  1994)  derived  a  control  input  as  a  function 
of  the  streamwise  wall  shear  only  by  modeling  the  near-wall  flow  with  a  spanwise 
velocity  growing  linearly,  and  normal  velocity  growing  quadratically,  with  normal 
distance  from  the  wall.  About  15%  drag  reduction  was  obtained  from  this  study. 
We,  however,  derive  a  simple  control  scheme  by  minimizing  cost  functionals  that 
are  related  to  the  streamwise  vortices,  which  have  been  found  to  be  responsible  for 
large  local  drag  in  turbulent  boundary  layers.  We  also  tried  to  minimize  drag  directly 
by  having  drag  itself  in  the  cost  functional,  but  it  was  not  successful  (see  §  3).  The 
objective  of  this  paper  is  to  demonstrate  that  a  wise  choice  of  the  cost  functional 
coupled  with  a  variation  of  the  formulation  can  lead  to  a  more  practical  control  law. 

We  present  how  to  choose  a  cost  functional  and  how  to  minimize  it  to  yield  simple 
feedback  control  laws  that  require  quantities  measurable  only  at  the  wall.  One  of 
the  laws  requires  spatial  information  on  the  wall  pressure  over  the  entire  wall  and 
the  other  requires  information,  also  over  the  entire  wall,  on  one  component  of  the 
wall  shear  stress.  We  then  derive  more  practical  control  schemes  that  only  require 
local  wall  pressure  or  local  surface  shear  stress  information,  and  show  that  they  work 
equally  well. 

2.  Suboptimal  procedure 

We  follow  a  similar  procedure  used  by  Choi  et  al.  (1993)  and  Bewley  &  Moin 
(1994).  The  problem  under  consideration  is  a  turbulent  channel  flow,  for  which  the 
governing  equations  are  the  Navier-Stokes  and  continuity  equations  with  the  no-slip 
boundary  condition: 
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whc:~  r  is  the  time,  xi.x2,x3  are  the  streamwise,  wall-normal,  and  spanwise  directions 
respectively,  u,-  are  the  corresponding  velocity  components,  p  is  the  pressure,  and 
Re  is  the  Reynolds  number,  and  the  control  input  is  the  wall-normal  velocity  at  the 
wall,  4 >- 

All  variables  are  non-dimensionalized  by  the  wall  shear  velocity,  uT,  and  the  channel 
half-width,  <5.  We  also  use  interchangeably  x,y,z  for  x,  and  u,v9w  for  uh  Periodic 
boundary  conditions  are  imposed  in  the  streamwise  and  spanwise  directions.  The  flow 
rate  in  the  streamwise  direction  is  kept  constant,  and  the  drag  is  measured  by  the 
mean  pressure  gradient  necessary  to  maintain  the  constant  flow  rate. 

We  found  that  the  choice  of  the  cost  functional  to  be  minimized  is  critical  in 
the  performance  of  the  control.  Since  the  streamwise  vortices  have  been  known  to 
be  responsible  for  large  drag  in  turbulent  boundary  layers,  we  tried  to  choose  the 
cost  functional  that  is  directly  related  to  them.  This  is  based  on  our  conjecture  that 
a  suitable  manipulation  of  the  streamwise  vortices  would  lead  to  drag  reduction. 
We  carefully  selected  two  cost  functionals  based  on  our  observation  of  a  successful 
control  of  Choi  et  al  (1994).  As  shown  in  figure  1,  Choi  et  al.' s  (1994)  blowing 
and  suction,  which  are  equal  and  opposite  to  the  wall-normal  velocity  component 
at  y+  =  10,  effectively  suppress  a  streamwise  vortex  by  counteracting  up-and-down 
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Figure  1.  Schematic  of  a  pressure  field  induced  by  a  control  based  on  y+  =  10. 
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Figure  2.  Correlation  between  a  pressure  field  of  a  no-control  case  and  a  pressure  field  modified 
by  a  control  based  on  information  at  y+  =  10.  A  line  indicating  no  change  of  the  pressure  field  is 
drawn  for  guidance. 


motion  induced  by  the  vortex.  This  blowing  and  suction  creates  locally  high  pressure 
in  the  near-wall  region  marked  with  *+’,  and  low  pressure  in  the  region  marked  with 
’  in  figure  1.  A  crucial  aspect  of  the  present  analysis  is  the  observation  that  this 
blowing  and  suction  increases  the  pressure  gradient  in  the  spanwise  direction  under 
the  streamwise  vortex  near  the  wall.  To  demonstrate  this  behaviour,  we  examine 
computed  flow  fields.  Figure  2  shows  a  scatter  plot  between  two  pressure  gradient 
fields:  dp/dzll,  is  the  pressure  gradient  before  the  control  is  applied  and  dp/dz\%  is 
the  pressure  gradient  at  the  same  location  after  the  control  of  Choi  et  al  (1994) 
is  applied  for  one  time  step.  It  is  apparent  that  the  control  increased  the  pressure 
gradient  significantly.  For  the  uncontrolled  case,  Re-  based  on  the  wall-shear  velocity 
and  the  channel  half-width  is  about  110.  The  spectral  code  of  Kim,  Moin  &  Moser 
(1987)  is  used  for  all  the  computations  presented  here.  Simulations  are  carried  out 
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using  32  x  65  x  32  spectral  modes  (with  dealiasing  in  the  streamwise  and  spanwise 
directions)  with  the  computational  domain  of  4nd  x  23  x  47r<5/3,  respectively. 

The  above  argument  suggests  that  we  should  seek  blowing  and  suction  that  increases 
the  pressure  gradient  in  the  spanwise  direction  near  the  wall  for  a  short  time  period 
(i.e.  in  the  suboptimal  sense)  in  order  to  achieve  a  similar  drag  reduction  to  that 
achieved  by  Choi  et  aV s  (1994)  control.  The  cost  functional  f(<p)  to  be  minimized  is 
then 

w-ssj.r^-'Vsj.r (*)'.**[  <i4) 

where  the  integrations  are  over  the  wall  (S)  in  space  and  over  a  short  duration  in  time 
At,  which  typically  corresponds  to  the  time  step  used  in  the  numerical  computation, 
and  £  is  the  relative  price  of  the  control  since  the  first  term  on  the  right-hand  side 
represents  the  cost  of  the  actuation  0.  Note  that  there  is  a  minus  sign  in  front  of 
the  second  term  since  we  want  to  maximize  the  pressure  gradient.  It  should  be  noted 
that  the  spanwise  pressure  gradient  at  the  wall  will  be  eventually  reduced  when  the 
strength  of  the  near-wall  streamwise  vortices  is  reduced  through  successful  control. 
Here,  blowing  and  suction  that  increase  the  spanwise  pressure  gradient  for  the  next 
step  are  sought  as  a  suboptimal  control. 

To  minimize  the  cost  functional,  we  first  define  the  differential  states  of  the  velocity 
and  pressure,  (0i?p),  using  a  Frechet  differential  (Finlayson  1972), 
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0  being  an  arbitrary  perturbation  field  to  0. 

'  Next,  we  choose  the  Crank-Nicolson  scheme  for  the  linea^  verms  and  a  Runge- 
Kutta  scheme  for  the  nonlinear  terms  to  yield  a  discretized  form  of  (2.1)  and  (2.2): 
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where  the  superscripts  n  + 1  and  n  denote  the  time  step,  and  Rn  includes  the  nonlinear 
terms  and  the  explicit  parts  of  the  pressure  gradient  and  viscous  terms.  The  Frechet 
differential  of  (2.8)-(2.10)  yields  the  governing  equations  for  the  differential  states 
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=  2-  (2-13) 
Note  that  (@Rn/@<j))  4>  =  0.  Hereinafter,  we  drop  the  superscript  n+ 1  and  all  variables 
are  understood  to  be  at  the  (n  +  l)th  time  step.  Note  that  there  is  no  contribution 
from  the  nonlinear  terms,  thus  making  the  equations  linear.  Generally,  the  suboptimal 
formulation  depends  on  the  time  advancement  scheme  used,  as  shown  here  (see  also 
Choi  et  al  1993).  Dropping  the  nonlinear  terms  may  miss  important  flow  dynamics. 
However,  we  found  from  our  numerical  tests  with  the  full  nonlinear  terms  included 
that  the  contribution  from  the  nonlinear  terms  is  negligible  in  our  boundary  control 
with  short  optimization  interval  At;  it  turns  out  that  the  conservation  of  mass  due  to 
the  wall  actuation  dominates  the  near-wall  dynamics. 

Under  the  approximation  that  2Re/At  >  fc2,  where  k  =  (fc2  +  k2)1/2,  and  kx  and 
kz  denote  the  streamwise  and  spanwise  wavenumbers  in  the  x-  and  z-directions 
respectively!,  the  above  equations  have  the  following  solutions  in  the  semi-infinite 
domain  with  periodic  conditions  in  the  x-  and  z-directions  (see  the  Appendix) : 
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where  9j,p,  and  $  are  the  Fourier  coefficients  of  8j,p>  and  f>  respectively.  For  the 
channel  geometry,  we  originally  considered  both  walls  and  found  that  the  interaction 
between  two  walls  is  negligible  as  long  as  the  typical  wavelength  (~  2n/k)  associated 
with  near-wall  structures  is  much  smaller  than  the  channel  width. 

The  Frechet  differential  of  the  cost  functional  (2.4)  becomes 
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The  Fourier  representation  of  the  above  equation  is 
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where  the  hat  denotes  the  Fourier  coefficient,  and  the  superscript  *  denotes  the 
complex  conjugate.  From  (2.17),  dp/dz  |w  can  be  expressed  in  terms  of  0\ 
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Equation  (2.19)  then  reduces  to 
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t  It  can  be  shown  that  2Re/Atk;iax  ~  Rel/ 4  >  1.  Here  we  used  uAt/Ax  —  0(1)  and 
kmax  ~  kn  ~  Re3/*7  where  kn  is  the  wavenumber  corresponding  to  the  Kolmogorov  length  scale. 
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Figure  3.  Correlation  between  a  spanwise  wall-shear  stress  of  a  no-control  case  and  a  spanwise 
wall-shear  stress  modified  by  a  control  based  on  y+  =  10.  A  line  indicating  no  change  of  the 
spanwise  wall-shear  stress  field  is  drawn  for  guidance. 


For  an  arbitrary  </>,  equation  (2.21)  should  be  satisfied,  yielding 
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From  the  requirement  that  the  Frechet  differential  of  the  cost  functional  be  minimized, 
i.e.  @Jf  /@4>  =  0,  the  optimum  4>  then  becomes 


4>  =  c|pw,  (2.23) 

where  C  is  a  positive  scale  factor  that  determines  the  cost  of  the  actuation.  Equation 
(2.23)  indicates  that  the  optimum  wall  actuation  is  negatively  proportional  to  the  sec¬ 
ond  spanwise  derivative  of  the  wall  pressure,  with  the  high-wavenumber  components 
reduced  by  1/k. 

Another  wall  quantity  that  indicates  similar  changes  of  the  near-wall  dynamics  due 
to  the  altered  pressure  field  is  the  spanwise  shear  at  the  wall,  dw/dy.  Owing  to  the 
added  pressure  gradient  in  the  spanwise  direction  below  the  streamwise  vortex,  the 
spanwise  flow  near  the  wall  is  also  induced,  thus  increasing  the  spanwise  shear  stress 
at  the  wall  (see  figure  1).  The  scatter  plot  between  the  spanwise  shear  stress  at  the 
wall  from  two  different  fields,  one  of  which  is  from  an  unperturbed  channel  and  the 
other  with  the  control  based  on  y+  =  10,  is  shown  in  figure  3.  Thus  another  choice 
for  the  cost  functional  to  be  minimized  is 
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Figure  4.  Correlation  between  a  pressure  field  of  a  no-control  case  and  a  pressure  field  modified 
by  a  control  based  on  equation  (2.23). 


Following  the  procedure  that  led  to  equation  (2.23)  yields  the  optimum  actuation: 
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Equation  (2.25)  indicates  that  the  optimum  wall  actuation  should  be  proportional  to 
the  spanwise  derivative  of  the  spanwise  shear  at  the  wall,  with  the  high-wavenumber 
components  reduced  by  1  fk.  Note  that  the  scale  factor  C  in  equations  (2.23)  and 
(2.25)  is  arbitrary.  As  mentioned  before,  the  suboptimal  formulation  depends  on  the 
time-advancement  scheme  used.  However,  if  we  used  a  different  implicit  scheme,  only 
the  resulting  constant  C  would  be  different.  This  does  not  cause  a  problem  since  C 
is  chosen  such  that  the  r.m.s.  value  of  the  wall  actuation  is  maintained  at  a  given 
value. 

In  the  following  simulations,  we  set  the  r.m.s.  value  of  (j>  to  be  equal  to  that  of  the 
wall-normal  velocity  at  y+  =  10,  which  gives  the  same  r.m.s.  value  of  wall  actuations 
as  that  of  Choi  et  al.  (1994). 


3.  Results 

The  above  control  laws  (2.23)  and  (2.25)  are  tested  in  a  turbulent  channel  flow. 
The  pressure  gradient  at  the  wall  is  monitored  to  see  if  the  equation  (2.23)  control 
increases  the  pressure  gradient  when  the  control  law  is  applied.  It  behaves  as  expected, 
as  shown  in  figure  4.  The  spanwise  shear  stress  at  the  wall  also  increases  when  the 
equation  (2.25)  control  is  applied  (see  figure  5).  These  results  also  confirm  that  a 
suboptimal  procedure  without  the  nonlinear  terms  does  not  cause  an  error  for  our 
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Figure  5.  Correlation  between  a  spanwise  wall-shear  stress  of  a  no-control  case  and  a  spanwise 
wall-shear  stress  modified  by  a  control  based  on  equation  (2.25). 


boundary  controls.  Since  these  control  laws  specifically  increase  the  wall  pressure 
gradient  or  spanwise  wall-shear  stress,  correlations  between  these  variables  before 
and  after  control  are  much  higher  than  those  corresponding  to  the  control  based  on 
y+  —  10  information.  This  suggests  that  the  controls  based  on  equations  (2.23)  and 
(2.25)  modify  the  flow  in  a  different  manner  from  the  control  based  on  information 
at  —  10. 

Time  histories  of  the  mean  streamwise  wall-shear  stress  f.r  different  control  laws 
are  shown  in  figures  6  and  7.  As  control  begins,  an  immediate  drop  in  the  shear  stress 
is  observed  for  all  cases.  At  the  same  expense  (i.e.  the  same  r.m.s.  value  of  $),  the 
control  based  on  (2.25),  which  reduces  drag  by  as  much  as  22%,  is  apparently  more 
effective  than  that  based  on  (2.23)  which  reduces  drag  by  16%.  This  indicates  that 
the  spanwise  wall-shear  stress  is  a  better  quantity  for  control  input. 

The  control  laws  presented  above,  however,  are  still  impractical  to  implement,  since 
they  are  expressed  in  terms  of  the  Fourier  coefficients  (i.e.  in  wavenumber  space), 
which  require  information  over  the  entire  spatial  domain.  Therefore,  the  inverse 
transforms  of  k;/k  and  i kjk  are  sought  numerically  so  that  the  convolution  integral 
can  be  used  to  express  the  control  laws  in  physical  space.  The  discrete  representation 
of  each  control  law  then  becomesf: 

<j)(XjyZk)  =  cEE^^xW’^'’  (11) 

/  k' 

t  Since  k:/k  and  i ks/k  are  not  periodic  in  the  wavenumber  space,  a  high  resolution  was  used  to 
obtain  Wpjk  and  WJ[  to  minimize  aliasing  error. 
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k(—  z/Az) 

0  1  2  3  4  5  6 

0  1.0000  —0.4427  0.0031  -0.0440  0.0044  -0.0138  0.0032 

j(=x/ Ax)  1  0.0413  -0.0007  -0.0040  -0.0037  -0.0029  -0.0021  -0.0016 

2  -0.0110  0.0057  0.0019  0.0011  0.0002  0.0000  -0.0003 

Table  1.  The  weight  distribution  of  equation  (3.1),  Wpjk.  The  weights  are  symmetric  between  j  and 

-;  and  k  and  -/c.  The  weights  are  normalized  by  W^.  Bold  faced  weights  are  used  in  the  calculation 
of  figure  6.  Here,  Ax+  =  40  and  A z+  =  13  (Ax  =  3Az). 


k{=  z/Az) 

0  1  2  3  4  5  6 

0  0.0000  1.0000  -0.1039  0.2679  -0.0852  0.1419  -0.0671 

j(=z  x/Ax)  1  0.0086  0.0537  0.0503  0.0310  0.0340  0.0148  0.0237 

2  0.0001  -0.0104  0.0059  0.0051  0.0100  0.0074  0.0092 

Table  2.  The  weight  distribution  of  equation  (3.2),  W?k.  The  weights  are  symmetric  between  j  and 
and  antisymmetric  between  k  and  —  k.  The  weights  are  normalized  by  W£v  Bold  faced  weights 
are  used  in  the  calculation  of  figure  7.  Here,  Ax+  =  40  and  Az+  =  13  (Ax  =  3Az). 


and 


dy 


(Xj+j'i  zk+k')t 


(3.2) 


respectively.  The  subscripts,  j  and  k,  denote  the  discretizing  indices  in  the  x-  and 
z-directions  respectively.  The  weights,  W?k  and  Wj are  given  in  tables  1  and  2.  Note 
that  the  weights  decay  rapidly  with  distance  from  the  point  of  interest,  suggesting  that 
the  optimum  actuation  can  be  obtained  by  a  local  weighted  average  of  the  pressure  or 
spanwise  shear  stress.  The  results  obtained  using  a  37-point  average  for  the  pressure 
and  an  11 -point  average  for  the  spanwise  shear  stress  yield  about  the  same  drag 
reduction  as  that  obtained  from  full  integration  using  equations  (2.23)  and  (2.25)  (see 
figures  6  and  7).  The  weights  used  in  these  calculations  are  bold-faced  in  tables  1 
and  2.  It  is  remarkable  that  localized  information  can  produce  such  a  significant  drag 
reduction,  especially  for  the  control  with  the  spanwise  shear  stress. 

The  weight  distribution  WJk  for  j  =  0  is  very  similar  to  the  one  found  in  the 
application  of  neural  networks  to  the  same  turbulent  flow  (Lee  et  al  1997),  in  which 
only  a  single  strip  of  the  spanwise  shear  information  was  used  in  the  network.  The 
control  law  of  Lee  et  al 


4>  = 


i/c,  dw 

K\  ? y 


w 


(3.3) 


produces  almost  the  same  distribution  of  blowing  and  suction  as  the  present 
results.  Note  that  blowing  and  suction  from  equation  (2.25)  are  nearly  the  same 
as  those  from  equation  (3.3)  because  the  near-wall  structures  have  relatively  slow 
variation  in  the  streamwise  direction  (the  kx  —  0  component  is  dominant).  This  blow¬ 
ing  and  suction  distribution  is  also  very  similar  to  the  one  based  on  y+  =  10  data 
(see  Lee  et  al.  1997). 

In  figure  8,  contours  of  the  streamwise  vorticity  in  a  cross-flow  plane  for  the  control 
with  equation  (2.25)  are  compared  with  a  no-control  case.  Significant  reduction  in 
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Figure  6.  Mean  streamwise  wall-shear  stress  for  various  control  laws  based  on  the  wall  pressure, 
compared  to  the  no-control  case:  thick  solid  line,  the  control  law  expressed  in  Fourier  space, 
equation  (2.23);  thin  solid  line,  the  control  law  expressed  in  physical  space,  equation  (3.1)  with  the 
integration  radius  of  6A2.  Ar+  ~  0.2. 


Figure  7.  Mean  streamwise  wall-shear  stress  for  various  control  laws  based  on  the  spanwise 
wall-shear  stress  compared  to  the  no-control  case:  thick  solid  line,  the  control  law  expressed  in 
Fourier  space,  equation  (2.25);  thin  solid  line,  the  control  law  expressed  in  physical  space,  equation 
(3.2)  with  1 1  points  in  the  spanwise  direction  only.  Ar+  ~  0.2. 


the  strength  of  the  streamwise  vortices  is  evident.  This  again  confirms  the  notion 
that  a  successful  manipulation  of  the  near-wall  streamwise  vortices  can  lead  to  drag 
reduction.  The  mean  streamwise  velocity  and  root-mean-square  velocity  near  the  wall 
are  shown  in  figures  9  and  10,  respectively,  and  compared  with  the  no-control  case. 
The  trends  are  very  similar  to  Choi  et  al  (1994)  and  Lee  et  al  (1997). 

Finally,  we  mention  that  we  tried  to  find  a  feedback  control  scheme  minimizing  the 
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(a) 


Figure  8.  Contours  of  the  streamwise  vorticity  in  a  cross-flow  plane:  (a)  no  control:  ( b )  control 
based  on  wall-shear  stress  (equation  (2.25)).  The  contour  level  increment  is  the  same  for  both  figures. 
Negative  contours  are  dashed. 


Figure  9.  The  mean  streamwise  velocity  normalized  by  ux :  Thick  solid  line,  control  law  (2.25); 
dashed  line,  no  control.  For  the  controlled  case,  the  velocity  is  normalized  by  the  controlled  ut. 
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Figure  10.  The  root-mean-square  velocity  normalized  by  u,  of  the  uncontrolled  flow:  thick  solid 
line,  control  law  (2.25);  dashed  line,  no  control. 


following  cost  functional: 


/(</>)  = 


2y4Ar 


p  rt+At 

Jsl  ' 


( p~  dr  dS  -|~ 


2AAt 


irm: 


d t  dS, 


(3.4) 


with  m  —  1  or  2.  This  cost  functional  is  the  most  natural  choice  since  it  contains 
the  quantity  directly  related  to  drag.  We  followed  the  same  procedure  to  derive  the 
following  control  schemes: 

$  =  0,  (3.5) 


for  m  —  1,  and 


i kx  du 

CT  8y 


H’ 


(3.6) 


for  m  —  2.  Equation  (3.5)  obviously  does  not  reduce  drag  and  equation  (3.6)  increased 
drag  in  our  numerical  simulations.  This  failure  appears  to  be  due  to  the  neglect  of 
the  nonlinear  terms  in  our  formulation,  since  the  numerical  solution  of  the  optimum 
wall  actuations  with  the  full  nonlinear  terms  gave  different  results  (Bewley  &  Moin 
1994).  This  suggests  that  even  in  a  short  time  interval  the  nonlinear  terms  should 
be  included  in  the  suboptimal  formulation  when  drag  itself  is  chosen  as  the  cost 
functional.  Manipulation  of  the  streamwise  vortices  by  wall  actuation,  however,  can 
be  accomplished  through  a  linear  process  as  shown  in  the  previous  section.  It  appears 
that  having  a  term  that  is  directly  related  to  near-wall  streamwise  vortices  in  the  cost 
functional  indirectly  includes  the  nonlinear  effect.  This  is  probably  related  to  the  fact 
that  near-wall  streamwise  vortices  are  self-sustained  by  a  nonlinear  process  (Hamilton, 
Kim  &  Waleffe  1995).  If  the  nonlinear  terms  are  included,  of  course,  it  is  impossible 
to  derive  the  feedback  control  law  in  closed  form  without  approximation.  Hill  (1993), 
on  the  other  hand,  considered  the  nonlinear  term  by  modelling  the  near-wall  flow 
using  the  Taylor  series  expansion  and  presented  a  feedback  control  law  in  closed 
form,  which  resulted  in  about  15%  drag  reduction. 
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4.  Summary 

We  have  obtained  two  simple  feedback  control  laws  for  drag  reduction  by  applying 
a  suboptimal  control  procedure  to  a  turbulent  flow.  This  was  possible  since  we  selected 
two  cost  functionals  guided  by  the  successful  control  based  on  quantities  monitored 
at  y+  =  10.  The  present  control  laws  perform  very  well  resulting  in  substantial 
drag  reductions  when  applied  to  a  turbulent  channel  flow.  More  convenient  control 
schemes  requiring  only  local  information  were  also  derived,  and  were  shown  to  work 
equally  well.  They  require  quantities  measurable  only  at  the  surface  and  thus  should 
be  easier  to  implement  in  practice.  The  present  results  further  substantiate  the  notion 
that  a  successful  manipulation  of  the  near-wall  streamwise  vortices  is  the  key  to 
boundary-layer  control  for  drag  reduction. 
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Appendix.  Derivation  of  the  solution  (2.14)-(2.17) 

Equations  for  the  differential  states  (0„p),  (2.11),  (2.12)  can  be  rewritten  in  terms 
of  the  Fourier  coefficients, 


At  (  d2 


2 Re  \  dy 


ti-£z  hb-*2  <U^p  =  o, 


ikxAt . 


'  2 Re  Vdy2  /  2  dy 


At  /  d2 


d  e, 

ikx6i  +  \k~62  H — ~ —  — 
dy 


with 


0,(0)  =  4>dn,  0/(  00)  =  0, 

where  0,  and  p  are  defined  as  follows: 

0,  =  ^]TeKy)e*‘*e*’% 


kx  k: 


p  =  Y.Y.p^'k’Xs'k::- 


kx  fc: 


(Al) 
(A  2) 
(A3) 
(A  4) 

(A  5) 

(A  6) 
(A  7) 


The  operation  i kx-  (Al)  +d/dy(A2)  +i k:-  (A3)  together  with  equation  (A 4)  yields 


d2p 


dy 


r  i  2  A  A 

--k-p  =  0, 


which  has  a  non-growing  solution, 


p  =  p*e  ky, 


(A  8) 

(A  9) 
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with  a  wall  value,  which  will  be  determined  later.  Using  (A  9),  we  can  find  solutions 
to  equations  (A  1),  (A 2),  (A3), 

e,(y)  =  yi M»  (exp[-(fc2  +  2Re/At)^y]  -  e^) ,  (A  10) 

03(y)=  yife.-Pw  (exp[-(fc2  +  2Re/Ar)1/:y]  -  e~ky) ,  (A  11) 

02(y)  =  ($~  y kP«j  exPE — (^2  +  2Re/At)1,2y]  +  y^e-^'.  (A  12) 


With  these,  equation  (A  4)  reduces  to 


exp[-(fc2  +  2Re/M)U2y]  =  0.  (A  13) 


Since  2Re/At  »  k2,  equation  (A  13)  yields 


2  i 

^  =  aTI0- 


(A  14) 


With  this,  equations  (A  10),  (All),  (A  12),  (A 9)  become  equations  (2.14),  (2.15),  (2.16), 
(2.17),  respectively. 
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stabilization  of  infinitesimal  and  finite-amplitude 
disturbances  in  plane  Poiseuille  flow 

By  SANJAY  S.  JOSHIf,  JASON  L.  SPEYER  and  JOHN  KIM 

School  of  Engineering  and  Applied  Science,  University  of  California,  Los  Angeles 
Los  Angeles,  CA  90024,  USA 

(Received  14  February  1995  and  in  revised  form  15  July  1996) 

A  systems  theory  framework  is  presented  for  the  linear  stabilization  of  two-dimen¬ 
sional  laminar  plane  Poiseuille  flow.  The  governing  linearized  Navier-Stokes  equations 
are  converted  to  control-theoretic  models  using  a  numerical  discretization  scheme. 
Fluid  system  poles,  which  are  closely  related  to  Orr-Sommerfeld  eigenvalues,  and 
fluid  system  zeros  are  computed  using  the  control-theoretic  models.  It  is  shown 
that  the  location  of  system  zeros,  in  addition  to  the  well-studied  system  eigenvalues, 
are  important  in  linear  stability  control.  The  location  of  system  zeros  determines  the 
effect  of  feedback  control  on  both  stable  and  unstable  eigenvalues.  In  addition,  system 
zeros  can  be  used  to  determine  sensor  locations  that  lead  to  simple  feedback  control 
schemes.  Feedback  controllers  are  designed  that  make  a  new  fiuid-actuator-sensor- 
controller  system  linearly  stable.  Feedback  control  is  shown  to  be  robust  to  a  wide 
range  of  Reynolds  numbers.  The  systems  theory  concepts  of  modal  controllability  and 
observability  are  used  to  show  that  feedback  control  can  lead  to  short  periods  of  high- 
amplitude  transients  that  are  unseen  at  the  output.  These  transients  may  invalidate  the 
linear  model,  stimulate  nonlinear  effects,  and/or  form  a  path  of ‘bypass’  transition  in  a 
controlled  system.  Numerical  simulations  are  presented  to  validate  the  stabilization  of 
both  single-wavenumber  and  multiple-wavenumber  instabilities.  Finally,  it  is  shown 
that  a  controller  designed  upon  linear  theory  also  has  a  strong  stabilizing  effect  on 
two  dimensional  finite-amplitude  disturbances.  As  a  result,  secondary  instabilities  due 
to  Liiinitesimal  three-dimensional  disturbances  in  the  presence  of  a  finite-amplitude 
two-dimensional  disturbance  cease  to  exist. 


1.  Introduction 

A  basic  problem  in  fluid  dynamics  is  the  theoretical  understanding  of  how  insta¬ 
bility  in  laminar  shear  flow  leads  to  transition  to  turbulence.  Since  laminar  flow' 
is  preferred  in  many  applications,  the  suppression  of  fluid  instabilities  to  maintain 
laminar  flow  would  be  very  useful.  Towards  this  goal,  active  boundary  layer  control 
of  instability  has  been  proposed.  The  nonlinear  aspects  of  the  transition  process  are 
still  not  clearly  known.  It  has  been  shown  (Orszag  &  Patera  1983)  that  some  shear 
flows  may  sustain  two-dimensional  finite-amplitude  instabilities  that  cause  infinites¬ 
imal  three-dimensional  disturbances  to  be  highly  unstable.  This  two-dimensional 
primary  instability /three-dimensional  secondary  instability  process  may  be  one  non¬ 
linear  mechanism  that  leads  to  transition.  Conversely,  the  behaviour  of  infinitesimal 

t  Present  address:  Jet  Propulsion  Laboratory,  4800  Oak  Grove  Drive,  Pasadena,  CA  91109,  USA. 
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(linear)  disturbances  in  laminar  shear  flow  is  well  understood  (Orr  1907;  Sommerfeld 
1908;  Drazin  &  Reid  1981).  Therefore,  linear  analysis  provides  a  natural  starting 
point  to  begin  to  develop  active  control  schemes  that  may  eventually  lead  to  full 
transition  control.  It  is  the  purpose  of  this  paper  to  develop  feedback  controllers 
based  on  linear  theory  that  stabilize  two-dimensional  plane  Poiseuille  flow  to  in¬ 
finitesimal  disturbances.  In  addition,  it  will  be  shown  that  a  controller  designed 
upon  linear  theory  has  a  strong  stabilizing  effect  on  two-dimensional  finite-amplitude 
instabilities.  As  a  result,  secondary  three-dimensional  instabilities  as  described  in 
Orszag  &  Patera  (1983)  cease  to  exist  in  such  a  system. 

Most  prior  work  in  the  area  of  laminar  flow  linear  instability  suppression  has 
concentrated  on  the  wave  superposition  approach.  A  nice  survey  of  past  work  is 
given  in  Joslin,  Erlebacher  &  Hussaini  (1994).  The  basic  idea  is  that  boundary  layer 
instabilities  appear  as  a  combination  of  many  sinusoidally  growing  waves  of  certain 
frequencies,  phases,  and  amplitudes.  If  these  wave  parameters  are  known,  or  if  they 
can  be  determined,  a  second  wave  may  be  stimulated  in  the  flow  that  is  exactly  out 
of  phase  with  the  instability  wave.  In  this  way,  the  two  waves  may  destructively 
interfere  and  flow  may  be  stabilized.  Disturbance  of  the  flow  field  to  cause  a 
wave  motion  to  appear  has  been  experimentally  demonstrated  using  several  methods 
such  as  vibrating  ribbons  (Schubauer  &  Skramstad  1947),  vibrating  wires  (Milling 
1981),  and  heating  elements  (Nosenchuck  1982).  In  addition,  several  authors  have 
numerically  obtained  wave  superposition  results  based  on  Navier-Stokes  simulations 
(Beringen  1984;  Metcalfe  et  al  1986).  Most  of  the  results,  however,  report  incomplete 
destruction  of  instability  waves.  Joslin  et  al.  (1994)  explain  that  wave  cancellation 
is  very  sensitive  to  the  wave  parameters  and  postulate  that  incomplete  destruction 
reported  in  past  studies  was  due  to  improper  phase  or  amplitude  properties  of  the 
cancelling  wave. 

In  addition,  Bower,  Kegelman  &  Pal  (1987)  considered  the  Orr-Sommerfeld  equa¬ 
tion  in  designing  an  input  to  channel  flow  that  may  counteract  the  effects  of  a 
disturbance  that  excites  instabilities.  They  showed  that  if  an  oscillating  flat  pulse 
function  at  the  lower  boundary  of  a  channel  excited  inherent  instabilities,  a  second 
oscillating  flat  pulse  function  could  be  constructed  downstream  that  negated  excite¬ 
ment  of  the  instabilities.  Much  like  the  wave  superposition  approach,  their  aim  was 
not  to  stabilize  the  underlying  dynamics  of  the  problem,  but  rather  to  ‘cancel’  the 
effects  of  a  specific  disturbance. 

Unlike  past  work,  our  aim  is  to  use  systems  theory  to  construct  a  combination  fluid- 
actuator-sensor-controller  system  that  is  inherently  stable.  In  essence,  the  approach 
changes  the  philosophy  of  the  problem  from  thinking  about  how  inputs  can  mitigate 
an  inherently  unstable  system  to  thinking  about  how  sensors  and  actuators  can  be 
added  to  form  an  entirely  new  stable  system.  Recently,  mathematical  control  systems 
theory  has  begun  to  be  applied  to  fluid  systems  (Burns  &  Ou  1994;  Gunzburger, 
Hou  &  Suobodny  1992;  Choi,  Moin  &  Kim  1993).  A  control  theory  approach  for 
laminar  flow  linear  instability  suppression  will  be  shown  to  have  several  advantages 
over  the  traditional  wave  cancellation  approach.  It  will  eliminate  the  need  to  explicitly 
measure  phase  and  frequency  of  instabilities.  Also,  it  will  provide  a  framework  to 
select  sensor  locations  in  order  to  have  the  least  complex  controller.  Further,  feedback 
control  will  be  shown  to  be  extremely  robust  to  changing  Reynolds  numbers  given 
proper  sensor  location.  In  addition,  linear  feedback  controllers  will  be  shown  to  have 
a  strong  stabilizing  effect  on  two-dimensional  finite-amplitude  disturbances. 

This  paper  is  organized  as  follows.  In  §2,  we  formulate  the  linear  channel  flow 
problem  using  the  linearized  Navier-Stokes  equations.  Boundary  input  in  the  form 
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of  blowing/suction  and  boundary  output  in  terms  of  shear  are  represented  within  the 
equations.  In  §3,  we  introduce  the  concepts  of  zeros  and  poles  of  a  system,  as  well  as 
control-theoretic  models.  The  governing  partial-differential  equation  for  the  system 
is  converted  into  a  set  of  first-order,  ordinary  differential  equations  via  a  Galerkin 
method.  These  first-order,  ordinary  differential  equations  are  then  converted  into  a 
control-theoretic  model.  Section  4  describes  the  infinite-dimensional  nature  of  the 
channel  flow  system  and  how  it  affects  the  selection  of  actuation.  Section  5  describes 
the  numerical  model  and  the  verification  of  the  calculated  poles  and  zeros.  Section  6 
introduces  feedback  control  design.  It  is  shown  that  judicious  sensor  placement,  based 
on  zero  locations,  can  lead  to  simple  control  schemes.  Furthermore,  the  control  system 
is  extremely  robust  to  change  in  Reynolds  number.  Section  7  explores  unobservable 
modal  reinforcement  as  a  possible  path  of  ‘bypass’  transition  in  a  controlled  system. 
It  then  shows  how  a  particular  control  model,  called  the  modal-canonical  state 
space  model,  may  be  used  to  assess  observability  of  modal  reinforcement.  Section  8 
demonstrates  multiple-wavenumber  instability  control.  Section  9  demonstrates  that  a 
linear  controller  has  a  strong  stabilizing  effect  on  two-dimensional,  finite-amplitude 
instabilities.  As  a  result,  secondary  three-dimensional  instabilities  as  described  by 
Orszag  &  Patera  (1983)  cease  to  exist.  Finally,  §10  outlines  conclusions. 


2.  Problem  formulation 


2.1.  Dynamic  equations 

We  consider  two-dimensional  plane  Poiseuille  flow  between  two  parallel  stationary 
plates.  Let  the  channel  be  of  finite  length  and  finite  height,  with  the  centreline  at 
zero.  The  flow  in  the  channel  is  described  by  the  unsteady  nonlinear  incompressible 
Navier-Stokes  equations.  In  order  to  study  the  linear  stability  of  the  system,  we 
follow  the  standard  procedure.  Consider  small  perturbations  in  the  velocities  of 
u(x9y,t)  in  the  horizontal  direction,  v(x,y,  t)  in  the  vertical  direction,  and  p(x9y,t)  in 
the  pressure  field.  Let  the  primary  flow  be  represented  by  U(y)  with  Uc  being  the 
centreline  velocity.  The  linearized  incompressible  Navier-Stokes  equations  may  be 
formed  by  substituting  the  primary  flow  and  small  perturbations  into  the  nonlinear 
incompressible  Navier-Stokes  equations  and  disregarding  the  second-order  terms 
involving  the  perturbations, 


£%i£l  +  uiy)^M  +  iHMst,.,-.,)  -  -2 +  ±v%w>, 

dt  dx  dy  dx  Re 


cv(x,y,t)  +  [/(v)^(x>y>r) 


dt 

du(x9y9t) 


cx 


dx 

6v{x,y,  t) 
dy 


_5p(x±t)  +  J_^y^ 

dy  Re 


=  0, 


(2.1) 

(2.2) 

(2.3) 


where  the  flow  variables  are  non-dimensionalized  by  the  channel  half-height,  H,  and 
centreline  velocity,  Uc.  Re  is  the  Reynolds  number  defined  as  (UcH/v)  where  v  is  the 
kinematic  viscosity.  By  introducing  a  ‘stream  function’,  y>(x,y,t). 


a  oip(x,y,t) 

u(x,y,t)  = - — - 

dy 


(2.4) 


and 


v(x,y,t)  =  — 


ci/;(x,>-,t) 

8x 


(2.5) 
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(2.1)— (2.3)  may  be  combined  into  a  single  equation, 


c  d2xp  5  d2\ p 
dt  ex1  +  dt  8y 2 


d2xp  ,  c  d2w 

-V^-^T,W  + 


d2U(y)dy 
dy2  dx  + 


iv2^- 


(2.6) 


Assume  periodic  boundary  conditions  in  the  streamwise  (x)  direction.  For  channel 
flow,  with  rigid  plates  at  y  —  — 1  and  y  =  1,  the  no-slip  boundary  conditions  become 


V>(x,y  =  —l,t)  =  0,  (2.7) 

^(x,y  =  -l,r)  =  0,  (2.8) 

dy 

tp(x,y  —  l,t)  =  0,  (2.9) 

^(x,y=l,r)  =  0.  (2-10) 

oy 

With  an  initial  condition 

w{x,y,t  =  0)  =  g(x,y)  (2.11) 


the  boundary  value  problem  is  completely  formed.  Existence  and  uniqueness  of  solu¬ 
tions  for  the  linearized  Navier-Stokes  equations  have  been  studied  in  Ladyzhenskaya 
(1969),  Kreiss  &  Lorenz  (1989),  and  Temam  (1984).  Equations  (2.6)-(2.11)  represent 
the  starting  point  for  construction  of  a  feedback  control  system.  These  equations 
neither  include  any  control  terms  nor  do  they  describe  any  sensing  of  flow  field 
variables. 


2.2.  Boundary  input 

We  consider  the  case  of  blowing/suction  at  the  lower  wall  of  the  channel.  The  bound¬ 
ary  conditions  are  now  modified  from  before  to  include  boundary  input,  represented 
as  the  known  separable  function  q(t)w(x)f(y ), 


v >(x,y  =  -Ft)  =  q(t)w(x)f(y  =  -l), 

(2.12) 

-zr-(x,y  -  -1,1)  -  q(t)w(x)  =  0, 

oy  oy 

(2.13) 

o' 

II 

II 

(2.14) 

X 

II 

II 

T 

X^ 

CU 

Oj  X; 

^  ii 

II 

O 

(2.15) 

Note  that  these  conditions  constrain  the  function  f(y)  such  that 

II 

JL 

iK 

o 

(2-16) 

8f(y  =  -l)  0 

dy 

(2.17) 

f(y  =  i)  =  o. 

(2.18) 

s/(y  =  i)  _  0 

dy 

(2.19) 

Many  functions  may  be  equally  appropriate.  One  such  function  is 

7(y)  =  V  +  !>'3-y2 -!>•  +  !• 

(2.20) 
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In  order  to  relate  boundary  conditions  on  xp  to  blowing/suction  in  the  wall-normal 
direction,  we  use  (2.5)  to  relate  v(x9y9t)  and  \p(x9y9t).  Then  (2.12)  becomes 

v(x,y  =  -1,  t)  =  -q(t)^ff(y  =  -1).  (2.21) 

Note  that  v{x,y,t\  is  related  to  the  derivative  of  w(x).  The  homogeneous  equation 
(2.6)  and  the  inhomogeneous  boundary  condition  (2.12)  can  be  converted  into  an 
inhomogeneous  equation  with  homogeneous  boundary  conditions  by  introducing 

<t>(x,  y,  t)  =  V(x,  y,  t)  -  q(t)f(y)w(x).  (2.22) 


Then  by  substituting  into  (2.6),  we  obtain 

5  82<t>  8  d2<l>  d3(l>  8  d24>  d2U{y)  8<j>  1  8A<j>  l  82  82<t> 

8t  ox2  +  8t  dy2  ~  (y)dx 3  (y)dx  8y 2  +  d y2  8x  Re  8x*  Redx 2  dy2 

+ 1% 

The  boundary  conditions  in  terms  of  <j>  are  now 


(2.23) 


<t>(y  =  —1)  =  0,  (2.24) 

o,  (2.25) 

dy 

4>{y  -  1)  =  0,  (2.26) 

(2.27) 

dy 


2.3.  Boundary  output 

We  use  the  streamwise  component  of  shear  at  a  single  boundary  point,  z(xi9y  —  —1,  t), 
as  our  boundary  output,  which  is  given  by 


z{xhy  =  — 1,  t)  = 


du(xi9y  =  -l,f) 
dv 


By  expressing  u(xi9y  = 


-1,0  in  terms  of  the  stream  function  (2.4), 
d2\p(xhy  =  -1,0 


z(x„y  =  —1,  r)  == 


dv2 


(2.28) 


(2.29) 


and  by  observing  (2.22) 

gV(x„y  =  -  1,Q  =  g^(x,-,y  =  -l,t)+  a2f(y  =  -l) 
dy2  ov-  dy - 

(2.30) 


z(xi9y  =  -l,t)  = 


3.  Zeros,  eigenvalues,  and  control-theoretic  models 

Linear  stability  analysis  of  (2.6)  (Drazin  &  Reid  1981)  shows  that  the  system  is 
linearly  unstable  for  a  range  of  Reynolds  numbers.  The  goal  of  this  paper  is  to 
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stabilize  the  system  using  control  theory.  To  do  this,  we  first  transform  the  governing 
equations  into  special  control-theoretic  models. 


3.1.  System  transfer  function 

Wc  have  defined  a  single-input/single-output  (SISO)  system  in  the  sense  that  only 
one  scalar  function,  q(t),  defines  the  entire  input  and  one  scalar  function,  z(t),  defines 
the  output.  A  common  form  of  control  model  for  a  finite-dimensional  SISO  system 
is  the  transfer  function  model.  The  transfer  function,  H(s ),  is  defined  as  the  Laplace 
transform  of  the  output,  z(t),  divided  by  the  Laplace  transform  of  the  input,  q(t ), 
where  zero  initial  conditions  are  assumed.  Then 


a  JSf[g(r)]  =  Z(s) 

*[«(*)]  fiW 


(3-D 


For  finite-dimensional  systems,  Z(s)  and  Q(s)  take  the  form  of  polynomials  in  the 
complex  variable  s.  These  polynomials  may  be  factored  to  yield  an  equivalent 
representation. 


H(s)  ^ 

H,s)  -  m 


J 


m*-w 

y=i _ 


ii(*- ft) 

i=i 


(3.2) 


In  this  form,  p, . .  .p,  are  the  poles  of  the  system.  The  poles  of  any  system  are  dependent 
solely  on  the  physics  of  the  underlying  system,  independent  of  any  particular  input  or 
output.  Unstable  modes  of  the  system  appear  as  poles  whose  real  part  is  greater  than 
zero.  As  will  be  seen  in  later  sections,  the  poles  are  closely  related  to  the  eigenvalues 
of  the  Orr-Sommerfeld  equation.  The  values  Ci  ■••0  are  the  zeros  of  the  system.  They 
are  heavily  dependent  on  which  particular  inputs  and  outputs  are  used  on  the  system. 
As  will  be  seen  later,  the  position  of  these  zeros  will  dictate  sensor  locations  and  will 
reveal  the  effect  of  feedback  control  on  the  eigenvalues.  A  graphical  representation  of 
the  transfer  function  can  be  produced  by  plotting  the  poles  and  zeros  in  the  complex 
s-s^ace. 


3.2.  State-variable  model 

Much  of  modern  control  theory  is  based  on  the  state-variable  representation  of  a 
dynamic  system.  This  representation  relies  on  the  basic  fact  that  the  motion  of  any 
finite-dimensional  dynamic  system  may  be  expressed  as  a  set  of  first-order  ordinary 
differential  equations.  As  a  simple  example  of  a  state  variable  model  (Franklin,  Powell 
&  Emami-Naeini  1988),  Newton’s  law  for  a  constant  single  mass,  M,  moving  in  one 
dimension,  x,  under  a  force,  F(t),  is 


M 


d2x(t) 

dt2 


=  F(0. 


(3.3) 


If.  we  define  one  state  variable  as  the  position  xf  =  x(t)  and  the  other  state  variable 
as  the  velocity  x^  =  dx(f)/df,  (3.3)  can  be  written  as 


dx2  F(t) 
d  t  ~  M  ' 


(3-4) 


(3-5) 
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Furthermore,  the  first-order  linear  ordinary  differential  equations  can  be  expressed 
using  matrix  notation 


dxi 

*dT 

'  0  1  : 

*1 

'  0  ‘ 

d*2 

0  0 

_  *2  _ 

i 

1 

.  IT  - 

^  =  Ax  +  Bq . 
at 

If  we  take  output  as  position,  Xi, 


(3.6) 


(3.7) 

(3-8) 


or 

z  =  Cx  (3.9) 

The  matrix  A ,  and  the  vectors  B  and  C  are  called  the  state  space  matrices  of  the 
single-input,  single-output  (SISO)  system.  More  specifically,  the  A  matrix  is  referred 
to  as  the  dynamic  matrix  of  the  system.  It  can  be  shown  that  the  poles  of  the  system 
are  simply  the  eigenvalues  of  the  A  matrix.  In  the  more  general  case  of  multiple-input, 
multiple-output  (MIMO)  systems,  B  and  C  are  matrices.  For  generality,  we  may  add 
another  term,  Dq,  to  the  output  equation  to  account  for  systems  in  which  there  is 
direct  feedthrough  from  the  input  to  the  output.  Then, 

z  =  C-x  +  Dq  (3.10) 


In  the  case  of  no  direct  feedthrough,  D  =  0,  the  state  space  and  transfer  function 
models  are  related  as 


m-W)’C(sl-AT 


(3.11) 


3.3.  State-space  formulation  for  channel  system 

We  will  convert  our  problem  into  a  set  of  first-order  ordinary  differential  equations 
and  then  form  a  state-space  model  from  these  equations.  The  state-space  model  can 
then  be  represented  with  transfer  function  poles  and  zeros.  We  will  proceed  in  this 
way  for  two  reasons.  First,  our  system  lends  itself  to  decomposition  into  first-order 
ordinary-differential  equations  by  use  of  a  Galerkin  method.  More  importantly, 
however,  the  state-space  model  lends  itself  to  extremely  elegant  ways  to  control 
eigenvalues  of  a  system  in  well  prescribed  ways.  It  should  be  noted  that  unlike  the 
single-mass  example  given  above,  the  channel  system  requires  an  infinite  number  of 
ordinary  differential  equations  to  describe  its  motion.  This  is  known  as  an  infinite- 
dimensional  system.  As  a  result,  any  finite  number  of  ordinary-differential  equations 
used  in  a  state-space  model  will  not  completely  describe  the  system.  The  difficulties 
associated  with  such  a  system  are  taken  up  after  a  discussion  of  the  Galerkin  method 
used  to  obtain  the  ordinary-differential  equations. 


3.3.1.  First-order  system 

We  use  a  standard  Galerkin  procedure  to  convert  the  governing  partial  differ¬ 
ential  equation  (2.23)  into  a  system  of  ordinary-differential  equations.  Assume  an 
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approximate  solution  of  (2.23)  as 


TV  AJ 

<t>a(x,y,  t)  =  m(t)P„(x)rm(y). 


n=—N  m= 0 


By  using  Fourier  functions,  e u'3cX,  for  P„(x)  and  basis  functions  constructed  from 
Chebyshev  polynomials  for  rm(y)  (Joshi  1996),  we  obtain  a  first-order  system.  Define 
inner  products  in  the  streamwise  (x)  and  normal  (y)  directions  respectively: 

[gmocx^Unoox^  A  1  ^  =  ^  Oo  =  ^,  (3-13) 

L  J-L/2  L 


r r  (  \  r  (  w  —  f  ^ n{y)^m{y) 

[rm(y)^n{y)]y  —  J  ^  ^  _  yiy/2  ^ 


—l  \  J  / 

where  L  is  the  non-dimensional  length  of  the  finite-length  channel.  Applying  the 
orthogonality  of  the  basis  functions  in  x,  we  obtain  a  system  of  first-order  ordinary 
differential  equations: 

-  E  +E^r&-E  -  E 

m=0  m=0  m= 0  m=° 

M  1  M  1  iU 

+  ^  aofe 

m=0  m=0  m=0 

+  ^U|„,(t) +  q(t){S[k},l  =  —N...N,k  =  0...M,  (3.15) 


:  =  -  c ,P,(x)  ^^,Ffc(>')  (3.16) 


^,P,(x)  [t/(yMy),A(.v)]y-  ^pP/(x)  U(y)^^,A(y) 


^.p,(*)l  f^(yxr*w 

5x  J  df 


4.  [^.wol  »),%)!, 

Re  tfX4 


+2^-[^,P((x)l  [^,r.(y)l  +^-Mx),P,(x)]J^,A(y)l  .(3.17) 
Re  [  dx2  J  x  [  cy-  J  v  Re  [  oy  J  y 

In  this  system,  a0  is  the  fundamental  wavenumber  in  the  x-direction,  defined  as  2n/L, 
and  the  /?  coefficients  are  defined  in  terms  of  the  following  scalar  products,  where  the 
F(y)  are  the  basis  functions  in  y: 

fL -  [^^r1-r‘w]j.  J- 0-".  <318> 

C  =  [U(y)rm(y),rk(y)]y,  (3.19) 

ft  -  [^r.W.r.(r)|  .  (J  ®) 
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(3.21) 


02+  A 
Pmk 


U(y)5-^^,rk(y) 


dy 2 


3.3.2.  Transformation  to  state-space  form 
We  may  visualize  (3.15)  in  block  matrix  form  as 

dai~-n 


r  Ms 


0 


0 


0  M^s+t  0  0 


0 


0 


0 


dt 

dfl;=-N+l 

d  t 


L  o 

0 

0 

MN 

J 

da/=tf 

L  dr 

K-n 

0 

0  1 

a\=-N 

0 

K-n+1 

0 

0 

0 

0 

0 

.  0 

0 

0 

K»  . 

+  [  r,  r2  ] 


<?(0 

<4(0 

dt 


(3.22) 


where  di=p  denotes  a  column  vector  of  all  a\k  coefficients  whose  first  index,  /,  is  the 
value  p.  In  compact  notation, 


M~=Ka+[Yl  Y2  ] 
dt  1  J 


9(0 
<4(0 
L  dt 


Assuming  non-singularity  of  the  M  matrix,  we  may  invert  to  obtain 


da 


—  =M~'Ka+  [  M~lYi  M~'Y2 
at 


q(t) 

dq{t) 

dr 


(3.23) 


(3.24) 


where  a,  K,  Y{  and  Y2  are  all  complex.  By  expanding  in  terms  of  real  and  imaginary 
parts,  we  obtain 

-----  -f  i—r—  —  M  1 K  rCLr  4~  iM  1  KfdR  4-  i M  RQ[  —  Ari  '/C/fl; 

dr  dr 

-YM~lY\Rq(t)  +  iM~lYuq(t)  4-  M~l  Y2r  ^  ^  4-  i  M~1Y2i 
where  the  subscript  R  or  I  indicates  real  or  imaginary  parts.  Define 


(3.25) 


-  A 

P  = 


aR 

ai 

q(t) 


Then  the  state  space  system  is 


dp 
dt 

We  write  (3.27)  as 


m-'kr 

m~'y]R  ' 

M  1Y2r 

M-'K, 

/W-’Kr 

M~'Yu 

P  + 

m-'y2I 

0 

0 

0 

1 

dq(t) 

dt 


^ f=Ap  +  B ^  =Ap  +  Bu(t) 
dr  dr 


(3.26) 


(3.27) 


(3.28) 
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where  the  input  to  the  system  is  dq(t)/dt.  Note  since  q(t)  in  (2.12)  is  related  to 
physical  blowing/suction,  it  is  important  to  keep  in  mind  that  the  input  to  this  model 
is  the  derivative  of  q(t).  In  physical  terms,  any  input  derived  from  this  model  will 
be  in  terms  of  dq(t)/dt  and  therefore  must  be  integrated  before  it  can  be  applied  as 
physical  blowing/suction. 


3.3.3.  Measurement  equations 

By  observing  (2.30)  and  using  the  assumed  solution  form  (3.12) 

,  ,  ^  52tp(x,,y  =  -l,r)  o2(j>  ,  ,32/(y  =  -i) 

=  - df -  =  Jyi  +  ^ - 3? — w(Xi) 

=  E  E aUtMxd81^:  ~1}  +  —Mxi)  +  p(t,Xhy  =  -1) 


n——N  m~ 0 


8v2 


dy2 


(3.29) 


where  p(t,  x„  y  =  —  1),  the  residual  component  of  streamwise  shear  due  to  the  truncated 
expansion,  is  assumed  to  be  negligible  compared  to  the  two  other  terms.  The  second 
term  on  the  right-hand  side  of  (3.29)  is  made  up  of  the  known  input  terms  (2.12). 
Denote 


cy2 


(3.30) 


By  pulling  out  the  complex  time  coefficients  and  denoting  them  as  the  vector  a  as 
above  we  may  construct  a  complex  observation  matrix,  0 


z(xhy  =  -l,t)  =  Oa  +  Dq(t). 


(3.31) 


Finally,  by  creating  p  by  stacking  the  real  and  imaginary  parts  of  a ,  as  well  as  q(t ), 
we  may  construct  an  observation  equation  in  state-variable  form 


'  zR(Xi,y  =  -l,t) 

■  Or  -O,  D  ' 

_  Zf(xt,y  =  — 1,  r) 

o 

o 

50 

O 

(3.32) 


Since  we  may  measure  only  real  output  (shear),  we  are  left  with  only  the  top  half  of 
the  observation  matrix, 


ZR(Xi,y  =  -l,f)  =  [  Or  -O,  D]p. 

In  order  to  describe  the  system  in  traditional  control  terms,  define 


(3.33) 


Then, 


c  =  [  Or  -O,  D], 
z*(Xi,y  =  -It)  =  op 


(3.34) 

(3.35) 


3.3.4.  Initial  conditions 

We  have  not  accounted  for  the  initial  value  (2.11)  in  our  boundary  value  problem. 
We  may  account  for  the  initial  condition  by  assuming  it  can  be  written  as  a  series 
expansion  in  terms  of  Pn(x)  and  rm(y).  Then 

<j)(xy y, t  =  0)  =  t p(x, y,r  =  0)  —  s(t  =  0)/ (y)w(x)  =  g(x, y)  - s(t  =  0)/(y)w(x) 

N  '  M 

=  E  Efc""p»wr»^) 

n——N  m= 0 


(3.36) 
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where  b ^  are  assumed  known.  Then 

a(t  =  0)  =  b.  (337) 

Since  a  and  b  are  both  complex,  stack  the  real  and  imaginary  components  to  arrive 
at  an  initial  condition  on  p : 


P(t  =  0) 


&R 

bi 

q(t  =  0) 


(338) 


We  have  now  defined  the  system  completely  in  the  state-variable  form  of  (3.7),  (3.9) 
where  A  and  B  are  given  in  (3.28),  C  is  given  in  (3.35),  and  the  initial  condition  is 
(338).  Although  we  have  shown  the  state-space  formulation  using  a  single  input  and  a 
single  output,  we  may  formulate  the  multiple  input/multiple  output  case  similarly.  In 
that  case,  the  multiple  input  functions  are  represented  as  a  sum  of  separable  functions 
in  (2.12)  and  multiple  outputs  are  represented  as  a  vector  whose  components  are  of 
the  form  (335).  It  should  be  noted  that  in  the  multiple  input  case,  it  is  important  for 
physical  implementation  that  for  any  finite  stretch  in  the  x-direction,  only  one  of  the 
multiple  input  functions  is  non-zero. 


4.  Infinite-dimensional  nature  of  the  channel  problem  and  effect  on 
actuator  design 

In  a  system  described  by  partial-differential  equations,  such  as  the  channel  flow 
problem,  no  finite-dimensional  model  will  capture  all  the  dynamics  of  the  system. 
Such  a  system  is  called  infinite-dimensional.  By  examining  the  Galerkin  approximate 
solution  (3.12),  we  see  that  the  approximate  solution  can  converge  only  as  N,  M  — ►  oo. 
This  would  result  in  an  infinite  number  of  ordinary  differential  equations.  Each  n 
value  included  in  the  approximate  solution  is  referred  to  as  adding  an  additional 
wavenumber  in  the  dynamic  model.  Even  for  just  one  wavenumber,  n  =  N  (say), 
we  see  the  assumed  Galerkin  solution  requires  an  infinite  number  of  terms  in  the 
y-direction, 

M— oc 

4>a(x,y,t)  =  ]T  aNm(t)PN(x)rm(y).  (4.1) 

m=0 

Interestingly,  by  examining  the  structure  of  the  A  matrix  in  the  resulting  state-space 
model,  we  see  that  the  dynamics  associated  with  each  individual  wavenumber  are 
decoupled  from  the  others.  This  is  characterized  by  a  block  diagonal  form  of  the  A 
matrix,  (4.2).  As  a  result,  we  may  separate  the  problem  of  determining  system  poles 
into  a  set  of  smaller  problems  that  include  only  one  wavenumber  at  a  time: 


A  = 


[4(71=1)]  0 

0  [4(rt=2)]  0 

0  0 

oo  0  0 


cc 

0 

0 

[4(n=oc)]  - 


(4.2) 


Indeed,  the  eigenvalues  (poles)  of  the  entire  A  matrix  are  simply  the  eigenvalues 
(poles)  of  each  smaller  A  matrix  block.  Offsetting  this  computational  advantage, 
however,  is  the  fact  that  the  locations  of  the  zeros  are  not  decoupled  by  wavenumber. 
Therefore,  even  though  the  poles  of  the  system  can  be  calculated  separately  using 
only  one  wavenumber  in  a  particular  computation,  the  zeros  of  the  system  force  all 
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the  wavenumbers  in  the  model  to  be  considered  together.  This  can  easily  be  seen 
from  an  examination  of  the  transfer  functions  associated  with  just  two  wavenumbers. 
Consider  two  transfer  functions  obtained  by  using  two  distinct  wavenumbers  and  the 
same  input,  u(t ),  and  output,  z(t): 

„  N(n=l)(s)  Tj  ty"-  2)(s) 

(n=1)“D<„=»(s)’  ,n=2)  Din=2)(sY 

The  roots  of  N(s)  are  the  zeros  and  the  roots  of  D(s)  are  the  poles.  In  terms  of  input 
and  output, 

z  =  H(„=i)U,  z  =  H(„=i)U.  (4.4) 

Considered  together,  however, 


Z  =  [#(/,=  1)  +  ff(n=2)]  M 


_  £>(n=2)(^Wn=l)(^)  +  Djn=l)(s)N(n=2)(s) 

D(n=l)(s)Di„=2)(s) 


The  poles  of  the  new,  combined  transfer  function  are  clearly  those  of  each  model 
separately.  The  zeros,  however,  are  the  roots  of  a  completely  new  polynomial 
£>0.=2)(s)N(n=i)(s)  +  D(n=l)(s)N^2)(s),  which,  in  general,  has  little  to  do  with  either  of 
the  original  numerators. 

We  have  seen  that  an  infinite  expansion  is  needed  to  fully  model  the  channel  system. 
Since  in  practice  an  infinite  expansion  cannot  be  used  to  create  a  state-space  model,  we 
will  always  have  some  part  of  the  system  dynamics  that  is  unmodelled.  Furthermore, 
we  have  seen  that  the  A  (dynamic)  matrix  may  be  decoupled  by  wavenumber. 

Consider  a  partition  of  the  state-space  model  as 


z  =  [  Cm  Cu  ]  x,  (4.7) 

where  .he  subscripts  m  and  u  represent  the  modelled  and  unmodelled  parts  of  the 
system  In  the  channel  flow  problem,  the  unmodelled  part  is  meant  to  denote  only 
the  dynamics  of  wavenumbers  left  out  of  the  finite-dimensional  model.  Note  that 
both  Au  and  Am  are  of  infinite  dimension;  Au  because  of  the  infinite  number  of 
wavenumbers  left  out  of  the  reduced-order  model  and  Am  because  of  the  infinite 
number  of  expansion  functions  needed  in  >'  for  each  of  the  finite  number  of  modelled 
wavenumbers.  One  way  to  avoid  considering  unmodelled  wavenumber  dynamics  is 
to  ensure  that  the  control  input,  u(t),  has  no  effect  whatsoever  on  the  unmodelled 
wavenumber  dynamics.  In  terms  of  the  state-space  model  (4.6),  this  is  equivalent  to 
rendering  Bu  =  0.  This  is  known  as  making  the  unmodelled  wavenumber  dynamics 
uncontrollable.  By  examining  (3.27),  we  see  that  the  B  matrix  is  formed  from  terms 
of  S,j  in  (3.16).  If  w(x)  in  (3.16)  is  such  that  S/k  =  0,  then  those  components  of  the 
B  matrix  are  zero.  Equivalently,  we  must  ensure  that  the  projection  of  w(x)  onto  the 
unmodelled  wavenumbers  is  zero.  Owing  to  the  orthogonality  of  Fourier  components, 
we  select  w(x)  to  be  made  up  of  modelled  Fourier  components  only. 


5.  Single-wavenumber  channel  model 

Consider  the  channel  model  shown  in  figure  1.  The  total  non-dimensional  length  of 
the  channel,  L,  is  47i.  In  this  case,  the  fundamental  wavenumber,  a0  =  0.5.  Only  one 
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M -  1  =  An  - ► 


Figure  1.  System  model  for  Poiseuille  flow.  Input  is  applied  along  the  bottom  plate  as 
w(x)  —  sin(x)  and  shear  is  measured  at  x  =  n/2  with  Re  =  10000  and  a  —  1.0. 


wavenumber  is  included  in  the  model,  corresponding  to  a  —  /cco  —  1.00.  The  Reynolds 
number  is  chosen  as  Re  =  10000.  Input  is  distributed  along  the  bottom  plate 
with  a  sinusoidal  weighting  function  in  order  to  render  the  unmodelled  wavenumber 
dynamics  uncontrollable  as  described  in  §4.  In  terms  of  the  input  function,  w(x),  see 
(212) 

w(x)  =  sin(x).  (5.1) 

Note  that  the  physical  blowing/suction,  T(x,y,£),  takes  the  form 


v{x,y,t)  =  -q(t  =  -1)  =  -q{t)  cos{x) f{y  =  -1)  (5.2) 

see  (2.21).  This  type  of  input  may  be  achieved  in  practice  by  a  large  number  of 
independently  controllable  actuators  that  are  distributed  along  the  lower  channel 
wall.  The  f{y)  function  is  chosen  as  in  (2.20).  The  sensor  location  is  x  =  +\n.  In 
order  to  visualize  the  control-theoretic  model,  the  system  A ,  B,  and  C  matrices  are 
transformed  to  transfer  function  form,  H{s).  Figure  2  shows  the  locations  of  the  poles 
and  zeros  in  the  5-plane  for  Me  single-wavenumber  model. 


5.1.  Relation  of  poles  to  eigenvalues  of  the  Orr-Sommerfeld  equation 

Poles  of  the  transfer  function  or,  equivalently,  eigenvalues  of  the  A  matrix  are  closely 
related  to  the  eigenvalues  of  the  Orr-Sommerfeld  equation.  Assume  a  solution  of 
(2.6)  of  the  form 

y>(x,y,t)  =  P(y)tixxc-'*a.  (5.3) 

By  substituting  into  (2.6),  we  obtain  the  familiar  Orr-Sommerfeld  equation  in  the 
normal-mode  form, 


[U(y)-c] 


my) 

8f- 


oy - 


l 

ictRe 


my) 

8y 4 


■  2a 


dy2 


+  “4 /?(>’)  • 
(5.4) 


Here,  the  Reynolds  number.  Re,  is  known;  the  wavenumber,  a,  is  assumed  real  and 
known;  the  complex  wave  speed,  c,  is  the  eigenvalue  of  the  problem;  and  the  function 
P(y)  is  the  eigenvector  of  the  problem.  Stability  of  a  flow  for  a  given  value  of  Re 
and  a  is  determined  by  the  imaginary  part  of  c.  If  the  imaginary  part  is  positive,  the 
solution  (5.3)  becomes  an  unbounded  exponential  and  flow  is  unstable.  Poles  of  the 
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Figure  2.  Pole  (x)  and  zero  (o)  configuration  for  the  system  model  of  figure  1.  Note  that  all  poles 
and  zeros  appear  in  complex-conjugate  pairs.  Each  pair  of  complex-conjugate  poles  represents  one 
mode  of  the  system.  One  pair  of  poles  is  to  the  right  of  the  imaginary  axis.  This  indicates  the  system 
has  one  unstable  mode.  Note  that  the  single  pole  at  the  origin  represents  a  built-in  integrator  due 
to  u{t)  -  dq{t)/dt.  Channel  model:  Re  =  10000,  shear  sensor  at  tt/2,w(x)  =  sin(x),  L  -  4n,  a  =  1.0. 


transfer  function  in  the  Laplace  domain  at  s  =  pt  correspond  to  solutions  in  the  time 
domain  of  ep>t.  As  a  result,  we  observe  that  poles  of  the  transfer  function  are  related 
to  eigenvalues  of  the  Orr-Sommerfeld  equation  as 

pt  =  — iac.  (5.5) 

In  order  to  validate  our  linear  code,  we  compare  our  poles  f  j  eigenvalues  produced  in 
Orszag  (1971).  Orszag  (1971)  obtained  Orr-Sommerfeld  eigenvalues  for  the  channel 
problem  with  a  =  1.00  and  Re  -  10  000.  He  reported  only  one  slightly  unstable 
eigenvalue  at  s  «  0.00373967  +  i0.23752649  (c  =  0.23752649  -  i0.00373967).  The 
eigenvalue  is  seen  as  unstable  by  its  positive  real  part.  We  obtained  identical  results 
including  the  one  unstable  mode  at 

5  =  0.00373967  +  i0.23752649. 

All  other  stable  modes  obtained  in  the  present  study  are  identical  to  those  reported 
in  Orszag  (1971).  The  goal  of  our  control  system  will  be  to  move  these  unstable  poles 
into  the  stable  half  of  the  s-plane  or,  equivalently,  make  sure  the  controlled-system 
poles  all  have  real  parts  less  than  zero. 

5.2.  Verification  of  model  zeros 

Verification  of  the  system  zeros  is  more  difficult  due  to  the  fact  that  no  published 
results  exist  to  our  knowledge.  We  use  the  channel  code  simulation  used  in  Kim, 
Moin  &  Moser  (1987)  to  verify  our  zeros.  This  code  is  a  spectral  channel  flow  code 
which  uses  periodic  boundary  conditions.  Consider  the  transfer  function  model  (3.2) 
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Figure  3.  Output  from  the  Navier-Stokes  simulation  with  input  applied  as  u(t)  =  e+0-31755^. 
Although  an  unbounded  input  is  applied,  near  zero  output  is  observed  at  sensor  location  n/2  (solid 
curve).  Sensor  location  0  (dashed  curve)  and  sensor  location  n  (dash-dot  curve)  do  not  show  zero 
output.  Channel  model:  Re  -  10000,  shear  sensor  at  it/ 2,w(x)  =  sin(x),  L  =  4rc,  a  =  1.0. 


with  J  =  1  and  1  =  1: 


In  the  Laplace  domain, 


H(s) 


Z(s )  (5 -Ci) 
U(s)  (s-piY 


sZ(s)  -  PiZ(s)  =  sU(s)  -  CiU(s). 


(5.6) 

(5.7) 


Assuming  zero  initial  conditions  and  transforming  into  the  time  domain,  (5.7)  becomes 


dz{t) 

“5 r 


-Piz(t)  = 


d  u(t) 
dt 


Ci«(r) 


(5.8) 


where  z(t)  is  the  output  function  and  u(r)  is  the  input  function.  If  the  input  is  taken 
as 

u(r)  =  ewr  .  (5.9) 

then  the  right-hand  side  of  (5.8)  becomes  zero  and  the  system  behaves  as  if  zero 
input  has  been  applied.  As  a  result,  z(t)  remains  zero  for  all  time.  Therefore,  we  may 
verify  zeros  of  the  system  by  applying  non-zero  input  in  special  ways  and  observing 
zero  output.  Note  from  figure  2  that  for  sensor  location  x  =  -F^rc,  one  zero  exists 
in  the  right  half  plane  at  5  =  +0.317557  +  Oi.  This  represents  an  ideal  zero  to  check 
as  this  corresponds  to  an  unbounded,  unstable  input.  Figure  3  shows  output  from 
the  Navier-Stokes  simulation  with  input  u{t)  =  dq(t)/dt  =  e+03I7557f.  Note  that  even 
though  an  unbounded,  exponentially  growing  input  is  applied  to  the  system,  near 
zero  output  is  observed  at  the  sensor,  thus  verifying  that  Ci  is  indeed  a  zero  of  the 
system.  This  is  contrasted  with  measurements  at  different  sensor  locations  that  grow 
rapidly.  With  a  sensor  at  a  different  location,  the  zeros  of  the  system  change  position 
so  zero  output  is  no  longer  expected  for  this  particular  input.  Absolute  zero  is  not 
observed  at  sensor  location  x  =  due  to  slight  numerical  inaccuracies  in  the  model. 


6.  Feedback  control  and  stabilization 

At  this  point,  with  the  construction  of  a  valid  state-variable  model,  any  number  of 
control  schemes  may  be  employed  to  stabilize  the  system.  A  general  control  system  is 
shown  in  figure  4.  The  output  of  the  system  is  fed  into  a  controller.  The  output  of  the 
controller  is  then  used  to  create  an  input  to  the  system.  The  design  of  a  controller  that 
achieves  certain  system  characteristics  is  the  goal  of  control  system  design.  Several 
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modern  control  techniques  may  be  applied  that  require  a  state-variable  model.  In. this 
paper,  the  primary  goal  will  be  system  stability.  For  this  purpose,  a  simple  constant 
gain  feedback  with  integral  compensator  will  be  shown  to  be  sufficient. 


6.1.  Constant  gain  feedback  with  integral  compensator 

Figure  5  shows  an  integral  compensator  feedback  control  scheme.  K  is  referred  to  as 
the  gain  of  the  feedback.  The  output  of  the  system,  in  this  case  shear,  is  multiplied 
by  a  feedback  gain,  integrated  in  time,  and  then  fed  back  as  blowing  and  suction  at 
the  input.  Note  that  the  integrator  is  required  because  of  the  structure  of  the  input 
given  in  (3.27),  i.e.  the  input  is  taken  as  the  derivative  of  suction/blowing.  The  signal, 
r(f),  is  called  the  reference  signal.  It  is  used  to  define  the  desired  output.  In  our  case, 
we  would  like  the  shear  output,  z(xhy  =  -1,0,  to  be  zero.  Therefore,  we  set  r(t)  to 
zero.  Then 

u(0  =  =  -fCz(x„j>  =  -1,0-  (6T) 

ot 

Therefore,  assuming  g(0)  =  0, 

<7(0  =  -K  [  z(Xi,y  =  -l,T)dT  (6.2) 

Jo 

and  by  observing  (2.21),  we  may  describe  physical  blowing  and  suction  at  the  bound¬ 
ary: 

v(x,y  -  -1,0  -  -q{t)~-^f{y  =  ”1)  =  Kcos(x)  [  z(xt-,y  =  -l,T)dr  (6.3) 

OX  Jo 


As  defined  earlier, 


Therefore, 


A  &[z(t)]  _  Z(s) 
<e[u{t)}  U(sY 

Z(s)  =  H(s)U(s). 


(6.4) 

(6.5) 


In  the  absence  of  feedback,  one  mode  is  unstable.  Also,  one  pole  exists  at  the  origin 
for  the  integrator.  Consider  a  new  transfer  function  from  the  reference  input,  r(f),  to 
the  output,  z(r),  in  the  presence  of  feedback: 


u(t)  =  r(t)  -  Kz(t). 


(6.6) 


By  taking  the  Laplace  transform  of  (6.6), 

U(s)  =  R(s)  -  KZ(s). 


Then 


Z(s)  =  tf(s)[K(s)-KZ(s)]. 


Finally, 

Z(s)  =__f/(s)_ 
R{s)  1  -KH(sY 
The  new  poles  of  the  feedback  system  are  defined  by 


1  —  KH(s)  —  0. 


(6.7) 

(6.8) 
(6.9) 

(6.10) 


As  K  gets  larger  and  larger,  it  is  clear  that  the  poles  of  the  new  system  tend  toward 
the  zeros  of  H(s).  In  this  way,  modes  of  the  system  can  sometimes  be  changed  to 
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Figure  4.  Feedback  control.  The  output  of  the  system  is  fed  into  a  controller  and  then  to  the  input. 


k  L  Integrator  Blowing/ 
suction 

* - 0* - 

Feedback  gain 

Figure  5.  Feedback  control  for  channel  system.  The  output  of  the  system  is  multiplied  by  a 
feedback  gain,  integrated  in  time,  and  then  fed  back  at  the  input.  The  reference  signal,  r(t), 
equals  0. 

form  a  stable  system  (Franklin  et  al  1988).  Unstable  modes  that  appear  as  poles 
on  the  right  hand  side  of  the  complex  s-plane  and  ‘marginally  stable’  poles  on  the 
Im(s)-axis  are  drawn  to  the  left-hand  side  by  applying  feedback. 

6.2.  Sensor  placement 

We  have  seen  that,  by  applying  feedback,  poles  of  the  system  will  eventually  be  drawn 
to  zeros  of  the  system.  In  the  channel  flow  system  of  figure  2,  any  feedback  will 
cause  either  the  pole  at  the  origin  (integrator  pole)  or  the  unstable  mode  poles  to  be 
drawn  to  the  zero  on  the  real  axis  in  the  right-hand  plane,  thus  making  the  system 
unstable  to  a  greater  degree.  Therefore,  finding  a  transfer  function  that  has  all  zeros 
in  the  left-hand  plane  becomes  an  important  objective.  The  poles  of  the  system  are 
independent  of  sensing  or  actuation.  However,  the  zeros  of  the  system  are  dependent 
on  both  the  type  and  location  of  sensing  and  the  type  and  location  of  actuation. 
Figure  6  shows  the  pole/zero  configuration  for  the  channel  model  with  the  shear 
sensor  at  three  different  locations.  Only  the  top  half  of  the  s-plane  is  shown  since 
the  bottom  half  is  a  mirror  image  projected  across  the  real  axis  as  shown  in  figure  2. 
The  poles  of  all  the  models  are  in  the  same  location  as  expected.  However,  the  zeros 
are  different  in  all  three  cases.  In  figures  6(a)  and  6(b) ,  we  observe  a  lone  zero  in 
the  right-hand  s-plane.  However,  when  the  sensor  is  placed  at  x  =  +7t,  figure  6(c) 
shows  all  zeros  in  the  left-hand  plane.  In  fact,  there  is  a  region  around  x  =  n  that 
results  in  all  zeros  in  the  left-hand  plane.  By  placing  a  shear  sensor  at  x  =  7t,  simple 
feedback  with  integral  compensation  will  allow  stabilization  with  the  proper  value 
of  gain,  K.  In  the  case  of  sensor  locations  that  result  in  right-hand-plane,  so  called 
‘non-minimum  phase’,  zeros  stabilization  is  still  possible.  However,  more  complex 
controllers  (Bryson  &  Ho  1975;  Ogata  1990)  must  be  designed  that  are  beyond  the 
scope  of  this  paper. 

6.3.  Root  locus  analysis  and  numerical  simulation 

One  way  to  visualize  how  system  poles  will  change  as  the  feedback  gain,  X,  changes 
is  to  construct  a  root  locus  plot.  This  is  a  plot  of  all  poles  of  a  system  as  the  feedback 
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Figure  6.  Pole  (x)/zero  (o)  configuration.  Channel  model:  Re  =  10000,  w(x)  —  sin(x),  L  —  4ti, 
a  =  1.0.  Only  the  top  half  of  the  s-plane  is  shown.  Shear  sensor  at  (o)  7t/4  ( b )  3n/4,  (c)  it. 


gain  varies  from  K  —  0  to  K  -  oo.  Figure  7  shows  such  a  plot  for  the  system  shown  in 
figure  6(c).  When  K  reaches  O  all  feedback  system  poles  (closed-loop-poles)  lie  on 
the  left-hand  plane  and  no  instabilities  exist  in  the  new  feedback  system.  Numerical 
results  obtained  from  the  Navier-Stokes  simulation  for  the  new  feedback  controlled 
system  are  shown  in  figure  8.  The  computation  is  carried  out  without  feedback  until 
r  =  50,  after  which  the  feedback  is  turned  on.  We  see  that  the  growing  instability  is 
quickly  suppressed.  At  the  instant  the  controller  is  turned  on,  the  simulation  shows 
a  high  transient  response  due  to  the  non-continuous  nature  of  the  input  at  that  time 
instant.  In  terms  of  the  state  space  defined  in  (3.28),  (3.35), 


^  =Ap  +  Bu(t), 

(6.11) 

z(xhy  =  — 1,  r)  =  Cp, 

(6.12) 

u(t)  =  —Kz(Xj,y  = 

(6.13) 

Then,  in  a  dosed  loop, 

^  =  Ap  -  B(Kz(xity  =  —l,t)) 

(6.14) 

=  Ap-  BKCp  =  (A-  BKC)p 

(6.15) 
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Figure  7.  Root  locus  plot  for  channel  system  with  input  as  in  the  channel  model  (figure  1)  and 
shear  output  at  x  =  n.  The  poles  start  at  the  open-loop  poles  shown  with  x’s.  A sK  increases,  they 
start  to  move  toward  the  location  of  the  system  zeros,  shown  as  o’s.  The  pole  at  position  s  =(0,0) 
moves  directly  to  the  left.  The  unstable  pole  moves  quickly  to  a  position  just  to  the  left  of  the 
imaginary  axis.  Near  0.7  <  y  <  1.0,  -0.6  <  x  <  -0.4,  we  see  a  pole  moving  towards  a  zero  that 
is  out  of  the  range  of  this  figure.  Channel  model:  Re  =  10000,  shear  sensor  at  n,  w(x)  =  sin(x), 
L  =  4n,  a  =  1.0.  Only  the  top  half  of  the  s-plane  is  shown. 


Figure  8.  Navier-Stokes  simulation  of  feedback  control.  Shear  is  multiplied  by  a  feedback  gain, 
integrated,  and  fed  back  into  the  input.  Channel  model:  Re  =  10000,  feedback  shear  sensor  at  n, 
w(x)  =  sin(x),  L  =  471,  a  =  1.0.  Curve  1  (solid)  shows  shear  at  location  n;  curve  2  (dot-dash)  shows 
shear  at  location  jr/2. 

The  solution  to  this  equation  is 

p(t)  =  tAt-BKCtp(t  =  0).  (6.16) 

This  solution  will  approach  zero  as  t  ->  a)  since  the  eigenvalues  of  A- BKC  (closed- 
loop  poles)  are  all  stable. 

6.4.  Robustness  in  the  presence  of  Reynolds  number  Uncertainty 
A  major  advantage  of  feedback  control  systems  is  their  robustness  to  system  un¬ 
certainty.  From  a  practical  point  of  view,  Reynolds  numbers  may  not  be  known 
exactly  or  may  change  frequently  as  in  the  flight  of  an  airplane,  for  example. 
Figure  9  shows  the  open-loop  pole/zero  configurations  for  systems  with  varying 
Reynolds  numbers.  Note  that  these  systems  start  with  all  zeros  in  the  left-hand 
s-plane.  A  root  locus  analysis  shows  that  a  feedback  system  with  K  =  0.1  stabilizes 
both  systems.  Indeed,  a  feedback  gain  of  K  =0.1  stabilizes  systems  for  a  wide 
range  of  Reynolds  numbers  from  Re  =  1000  to  Re  =  40  000.  Figure  10  shows 
the  least-stable  pole  in  both  the  controlled  and  uncontrolled  systems  for  several 
Reynolds  numbers.  Recall  that  a  pole  with  real  part  less  than  zero  is  stable.  We 
see  near  Re  =  5772,  an  unstable  pole  appears  in  the  open-loop  (uncontrolled) 
system.  Unstable  eigenvalues  continue  to  exist  in  the  open-loop  system  until 
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Figure  9.  Pole  (x)/zero  (o)  configuration.  Channel  model:  shear  sensor  at  n,  w(x)  =  sin(x), 
L  =  47:,  a  =  1.0.  Only  the  top  half  of  the  5-plane  is  shown.  ( a )  Re  =  7500,  (b)  Re  =  12  500. 


Figure  10.  Real  part  of  least-stable  pole  for  open-loop  (solid  line)  and  closed-loop  (dashed  line) 
channel  system  with  1000  ^  Re  ^  40000.  Shear  sensor  at  it,  w(x)  =  sin(x),  L  -  4tt,  a  =  1.0.  Data 
points  are  shown  as  x’s. 


30000  ^  Re  ^  35  000.  An  identical  feedback  controller,  with  gain  K  =  0.1, 
however,  stabilizes  the  system  for  Reynolds  numbers  in  the  range  1000  ^  Re  ^ 
40000.  Clearly,  the  feedback  controller  is  extremely  robust  to  changes  in  Reynolds 
number. 
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7.  Unobservable  transient  response  as  a  path  of  bypass  transition 

Recently,  some  authors  (Trefethen  et  al  1993;  Butler  &  Farrell  1992;  Farrell  1988; 
Henningson  1994)  have  suggested  a  possible  path  of  bypass  transition  in  Poiseuille 
flow  that  is  caused  by  high  transient  response  due  to  the  non-self-adjointness  of  the 
evolution  matrix.  This  high-transient  behaviour  is  due  to  the  non-orthogonality  of  the 
eigenmodes  and  not  caused  by  any  single  mode,  but  rather  a  combination  of  many 
modes.  In  this  section,  we  suggest  a  possible  different  form  of  bypass  transition  in  a 
controlled  system  due  to  a  single  mode.  In  the  presence  of  any  type  of  control,  it  is 
possible  that  energy  from  an  input  is  fed  into  a  mode  such  that  the  mode  is  reinforced 
during  the  transient  period  before  it  eventually  dies  out.  During  the  transient  period, 
however,  the  mode  may  grow  to  an  amplitude  large  enough  to  trigger  nonlinear 
effects  that  induce  transition  to  turbulence.  This  may  be  a  possible  path  of  bypass 
transition  in  a  controlled  system. 

Single  modes  may  be  characterized  in  linear  systems  in  terms  of  modal  controlla¬ 
bility  and  modal  observability.  Modal  controllability  implies  that  a  particular  mode 
may  be  affected  by  the  actuators  chosen  (blowing/suction  in  our  case).  Modal  ob¬ 
servability  implies  that  a  particular  mode  may  be  measured  by  the  sensors  chosen 
(shear  in  our  case).  A  particular  mode  may  be  non-observable,  non-controllable,  or 
both.  This  will  be  seen  in  §§7.1  and  7.2.  If  high-transient  modes  are  controllable  and 
observable,  control  theory  may  be  used  to  suppress  them.  However,  if  high-transient 
modes  are  unobservable,  the  high  transients  will  never  be  seen  and  feedback  control 
cannot  be  used  to  suppress  them. 

7.1.  Modal  canonical  form 

The  state-space  formulation  provides  an  excellent  framework  for  assessing  the  rein¬ 
forcement  of  each  mode  individually.  Consider  an  n- dimensional  state  space  with 
scalar  input  and  output 

^ =  Ax  +  Bu ,  (7.1) 

dr 

z  =  Cx.  (7.2) 

We  may  perform  a  similarity  transformation  on  the  system  to  produce  a  new  state- 
space  representation  with  the  same  input-output  relationship,  but  with  a  different 
interpretation  of  the  state  variables.  Let 

P=[v,  v2  ...  v„  ]  (7.3) 

where  vj  is  the  ;th  eigenvector  of  the  A  matrix  and  n  is  the  dimension  of  the  A  matrix. 
A  new  representation  is  constructed  as  (Grace  et  al  1992;  Kailath  1980) 

^  =p-lAPx  +  P~'Bu,  (7.4) 

dt 

z  =  CPx  (7.5) 

where  x  =  P~‘x.  The  new  representation  is  written  as 

^  =  Ax  +  Bu ,  (7.6) 

dt 

z  =  Cx  (7.7) 

where  the  real  eigenvalues  of  the  original  A  matrix  appear  on  the  diagonal  of  A  and 
the  complex  eigenvalues  appear  in  a  2  x  2  block  on  the  diagonal  of  A.  For  an  A 


178 


S.  S .  Joshi,  J.  L.  Speyer  and  J .  Kim 


matrix  with  eigenvalues  (olu<t  ±;cd,a2),  the  4  matrix  is 


ai  0  0  0“ 

0  a  co  0 

0  —co  a  0 

.  0  0  0  a2  _ 


(7.8) 


In  this  form,  each  state-variable  pair  represents  a  mode  of  the  system.  Furthermore, 
the  modes  are  block  decoupled  so  that  each  mode  is  represented  by  the  2  x  2  system 


dxi 

1 7 

dx2 

dt 


a 

—co 


z\  =  [  Ci 


CO 

<7 


Ct 


*1 

*2 


*1 

X? 


51 

52 


U, 


(7.9) 

(7.10) 


where  co,  cr,  S2,  Ci>  and  c2  are  all  scalars.  This  is  known  as  the  modal  canonical 
state-space  form.  All  modes  can  be  monitored  directly  for  high-transient  reinforcement 
in  this  form. 


7.2.  Modal  observability  and  controllability 

It  is  easy  to  see  that  if  b\  and  b2  are  zero,  no  input  can  affect  the  mode  and  the 
mode  is  uncontrollable.  Similarly,  if  cx  and  c2  are  zero,  then  the  motion  of  xx  and 
x2  cannot  be  measured  at  the  output,  zu  and  the  mode  is  unobservable.  In  terms 
of  stabilization,  an  uncontrollable  mode  is  not  problematic  as  long  as  it  is  stable. 
However,  if  an  unstable  mode  is  uncontrollable,  nothing  can  be  done  to  stabilize  it. 
Observability  may  be  a  problem  even  if  modes  are  stable.  If  an  unobservable  mode 
is  highly  amplified  by  control  energy,  no  attempt  could  be  made  to  suppress  it  since 
it  would  not  be  observed  at  the  output.  In  terms  of  poles  and  zeros,  uncontrollable 
or  unobservable  modes  both  show  up  as  pole/zero  cancellations  in  the  complex  s- 
planc.  As  can  be  seen  from  figure  2,  many  poles  and  zeros  for  the  channel  system 
lie  on  top  of  each  other,  indicating  that  certain  modes  in  the  system  are  either 
^.controllable,  unobservable,  or  both.  Although  a  mode  may  not  be  physically 
observable  at  the  output,  the  modal  canonical  formulation  allows  us  to  numerically 
observe  the  state  evolution  of  each  mode  directly.  In  this  way,  we  may  assess  the  risk 
of  highly  amplified  modes  triggering  bypass  transition.  This  is  done  for  the  controlled 
system  simulated  in  figure  8.  The  uncontrolled  linear  system  is  simulated  with  an 
initial  condition  for  200  time  steps  at  which  time  a  feedback  controller  with  gain 
K  —  0.1  is  turned  on.  Figure  11  shows  the  state  evolution  of  the  modes  with  poles 
at  s  =  -0.1474  ±  0.8514i  and  s  =  -0.3252  ±  0.6361i  as  well  as  the  state  evolution  of 
the  one  unstable  mode.t  In  addition,  figure  11(a)  shows  the  shear  measurement  at 
the  output.  We  see  that  the  unstable  mode  dies  out  quickly  as  soon  as  the  feedback 
is  turned  on.  In  addition,  the  mode  at  5  =  -0.1474  ±  0.8514i  gains  almost  no  energy 
after  the  feedback  is  activated.  This  indicates  that  the  mode  is  not  a  bypass  mode. 
The  mode  at  s  =  —0.3252  ±0.636Ii  is  seen  to  gain  a  lot  of  energy  after  the  controller 
is  started.  Indeed,  the  amplitude  of  the  transient  response  is  nearly  twice  that  of 
the  unstable  mode.  This  represents  a  possible  bypass  mode  since  such  an  amplitude 

t  Note  that  amplitudes  shown  as  a  result  of  linear  system  simulations  should  only  be  used  for 
comparison  with  each  other  since  the  linear  model  has  been  scaled  to  unity  open-loop  feedforward 
gain,  i.e.  H{s )  =  kZ{s)/U{s)  where  k  =  1. 
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Time  steps  (mode  s  =  -0.3252+0.636 li) 


Figure  11.  Linear  model  state  evolution  of  system  modes  with  control  applied  at  r  =  200:  (a) 
measured  shear,  (6)  the  state  evolution  of  the  unstable  mode,  (c)  the  state  evolution  of  the  mode 
at  5  =  —0.1474  +  0.8514i,  (d)  the  state  evolution  of  the  mode  at  5  =  — 0.3252  +  0.636  li.  Channel 
model:  Re  —  10000,  shear  sensor  at  rc,  w(x)  =  sin(x),  L  =  4?r,  a  =  1.0. 


may  force  the  system  into  the  nonlinear  region  where  transition  to  turbulence  may  be 
triggered.  Furthermore,  the  high  energy  nature  of  the  mode  cannot  be  observed  at  the 
shear  sensor  since  this  mode  is  unobservable.  This  can  be  seen  from  the  plot  of  the 
shear  sensor  output.  Modal  un-observability  is  also  suggested  by  noting  the  relatively 
low  ci  and  c2  values  of  0.5155,  -  -0.1131  compared  to  the  c\  and  c2  values  of  —162.85, 
234.66  for  the  unstable  mode  .which  is  clearly  seen  at  the  output.  Fortunately,  the  high 
energy  nature  of  the  mode  ^es  not  lead  to  a  bypass  transition  in  this  case.  This  is 
verified  by  the  Navier-Stokes  simulation  in  figure  8,  which  shows  no  nonlinear  effects. 
If  the  instability  were  allowed  to  grow  to  a  higher  amplitude  before  the  controller 
was  applied,  however,  the  highly  amplified  mode  in  the  transient  response  might  have 
triggered  nonlinear  effects. 


8.  Multiple  instability  control 

Although  most  work  has  focused  on  suppression  of  single  instabilities,  more  realistic 
models  should  include  multiple  wavenumbers.  Indeed,  for  a  given  Reynolds  number, 
an  infinite  number  of  wavenumbers  exist.  Each  wavenumber  contributes  its  own  poles 
and  zeros  to  the  control-theoretic  model.  Consider  a  model  with  non-dimensional 
channel  length  2071,  where  input  is  applied  as  w(x)  =  sin(x)  -F  sin(0.9x).  In  this  model, 
wavenumbers  of  0.9  and  1.0  are  included.  Both  wavenumbers  lead  to  unstable  modes. 
The  pole/zero  configuration  of  this  new  two-wavenumber  model  with  shear  sensor 
at  7i  is  shown  in  figure  12(a).  Note  that  in  comparison  to  the  one-wavenumber 
model,  more  poles  and  zeros  exist.  The  original  one-wavenumber  model  contained 
one  ‘fork1  structure  of  poles,  while  the  two-wavenumber  model  contains  two  ‘fork’ 
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Figure  12.  (a)  Pole  (x)  and  zero  (o)  configuration  of  channel  system  of  length  20tt,  including 
wavenumbers  of  0.90  and  1.00.  Re  =  10000,  shear  sensor  at  tt,  w(x)  =  sin(x)  +  sin(0.9x).  Only  the 
top  half  of  the  s-plane  is  shown.  ( b )  Same  as  (a)  but  showing  closed-loop  poles  after  feedback  with 
gain  K  =0.1. 


Figure  13.  Navier-Stokes  simulation  of  channel  system  of  length  207:,  including  wavenumbers  of 
0.90  and  1.00.  Re  =  10000,  feedback  shear  sensor  at  n,  w(x)  =  sin(x)  +  sin(0.9x),  feedback  gain 
K  =0.1.  Curve  1  (solid)  shows  shear  at  location  71;  curve  2  (dash-dot)  shows  shear  at  location  2 n. 


structures.  It  can  be  visualized  that  in  a  model  with  several  wavenumbers,  several 
‘forks’  will  stack  on  top  of  each  other  in  the  s-plane.  Near  s  =  i0.2  are  two 
unstable  poles  in  the  right-hand  plane  that  represent  the  two  unstable  modes  in  the 
system.  As  seen  before,  all  zeros  lie  in  the  left  hand  s-plane.  Figure  12(b)  shows 
the  closed-loop  poles  after  feedback  with  gain  K  =  0.1.  Results  from  the  Navier- 
Stokes  simulation  are  shown  in  figure  13.  The  computation  is  carried  out  without 
feedback  until  t  =  200.  A  combination  of  two  growing  waves  is  seen  at  the  shear 
sensor.  At  t  =  200,  feedback  with  integral  compensation  is  applied  and  the  system  is 
stabilized. 
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9.  Effect  of  linear  controllers  on  two-dimensional  finite-amplitude 
disturbances 

The  critical  Reynolds  number  above  which  linear  instabilities  exist  in  plane 
Poiseuille  flow  is  Re  «  5772.  Below  this  Reynolds  number,  Poiseuilie  flow  is  linearly 
stable.  However,  transition  from  laminar  to  turbulent  flow  occurs  in  experiments 
at  much  lower  Reynolds  numbers,  typically  around  Re  =  1000.  Moreover,  even 
for  super-critical  Reynolds  number  flows,  linear  instabilities  predict  extremely  slow 
growth  rates.  In  practice,  transition  occurs  orders  of  magnitude  more  quickly.  There¬ 
fore,  linear  stability  alone  does  not  dictate  transition.  Transition  to  turbulence  is 
generally  accepted  to  be  a  nonlinear  three-dimensional  phenomenon.  In  previous  sec¬ 
tions,  we  constructed  two-dimensional  linear  controllers  based  on  the  two-dimensional 
linearized  Navier-Stokes  equations.  By  applying  these  controllers  to  plane  Poiseuille 
flow  with  infinitesimal  two-dimensional  disturbances,  we  saw  that  we  may  linearly 
stabilize  the  system  so  that  the  flow  could  no  longer  support  those  disturbances.  If 
we  apply  our  linear  controllers  to  flows  that  contain  finite-amplitude  disturbances,  we 
can  no  longer  use  linear  analysis  to  describe  the  flow  dynamics  as  nonlinear  terms 
become  relevant.  However,  the  application  of  a  controller  based  on  linear  analysis 
does  change  the  nonlinear  system. 

Many  authors  (Orszag  &  Patera  1983;  Bayly,  Orszag  &  Herbert  1988)  have  explored 
the  effect  of  two-dimensional  finite-amplitude  disturbances  on  three-dimensional  in¬ 
finitesimal  disturbances  in  plane  Poiseuille  flow.  Above  Re  «  2900,  for  two  dimensions, 
stable  non-attenuating  finite-amplitude  equilibria  exist  in  plane  Poiseuille  flow.  Below 
Re  «  2900,  non-decaying,  finite-amplitude  equilibria  do  not  exist.  However,  in  flows 
with  1000  ^  Re  ^  2900,  the  timescale  for  decay  is  so  large  that  the  flow  may  be 
considered  in  ‘quasi-equilibrium’.  It  has  been  shown  that  in  the  presence  of  such 
two-dimensional  finite-amplitude  disturbances,  infinitesimal  three-dimensional  distur¬ 
bances  are  highly  unstable  and  may  cause  transition  in  shear  flows.  Figure  14  shows 
the  energy  of  a  single-wavenumber,  two-dimensional  finite-amplitude  disturbance 
(a  —  1.0)  and  a  single-wavenumber-pair,  three-dimensional  infinitesimal  disturbance 
(a  =  1.0,/?  =  +1.0)  obtained  through  direct  numerical  simulation  at  Re  =  3000.  The 
maximum  amplitude  of  the  two-dimensional  finite-amplitude  disturbance  is  0.1  Lc. 
At  this  Reynolds  number,  wavenumber,  and  initial  energy  level,  we  see  the  tri¬ 
dimensional  finite-amplitude  disturbance  decaying  slowly.  The  three-dimensional  dis¬ 
turbance,  on  the  other  hand,  is  seen  to  rapidly  gain  energy.  Orszag  &  Patera  (1983) 
show  that  the  two-dimensional  instability  acts  as  a  mediator  for  transfer  of  energy 
from  the  mean  flow  to  the  three-dimensional  disturbance,  but  does  not  directly  provide 
energy.  The  growth  rate  of  the  three-dimensional  disturbance  is  orders  of  magnitude 
larger  than  that  of  the  Orr-Sommerfeld  instabilities. 

Orszag  &  Patera  show,  in  the  case  where  Re  ^  1000,  that  the  attenuation  of 
finite-amplitude  two-dimensional  disturbances  is  large  enough  that  three-dimensional 
disturbances  do  not  become  unstable.  The  fact  that  high  attenuation  of  the  two- 
dimensional  finite-amplitude  disturbance  prevented  three-dimensional  instability  be¬ 
low  Re  «  1000  suggests  that  if  controllers  can  be  created  that  speed  up  the  attenuation 
of  the  finite-amplitude  two-dimensional  disturbance  for  flows  with  Reynolds  num¬ 
bers  greater  than  1000.  three-dimensional  instability  may  be  eliminated.  We  have 
seen  in  §6  that  linear  controllers  did  stabilize  infinitesimal  two-dimensional  distur¬ 
bances.  Figure  15  shows  the  effects  of  the  linear  controller  of  §6  when  applied 
to  a  system  with  the  same  finite-amplitude  two-dimensional  disturbance  and  in¬ 
finitesimal  three-dimensional  disturbance  shown  in  figure  14.  The  linear  controller 
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Figure  14.  Energy  of  a  two-dimensional  finite-amplitude  disturbance  (a  =  1.0)  and  a 
three-dimensional  infinitesimal  disturbance  (a  =  1.0,/?  —  4-1-0)  for  Re  =  3000.  Solid  line  rep¬ 
resents  the  total  energy  of  the  two-dimensional  finite-amplitude  disturbance  and  the  dotted  line 
represents  that  of  the  three-dimensional  infinitesimal  disturbance. 


Figure  15.  Effect  of  linear  controller  when  applied  at  t  =  15  for  Re  —  3000.  The  top  two  lines 
show  the  two-dimensional  finite-amplitude  disturbance  without  (solid)  and  with  (dash-dot)  linear 
control.  We  see  at  t  =  1 5  that  the  controlled  two-dimensional  finite-amplitude  disturbance  is  highly 
damped.  The  bottom  two  lines  show  the  infinitesimal  three-dimensional  disturbance  without  (dot) 
and  with  (dashed)  linear  control.  That  without  controller  gains  energy  quickly,  while  that  with 
controller  quickly  fails  to  gain  energy. 
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dramatically  increased  the  attenuation  of  the  finite-amplitude  two-dimensional  dis¬ 
turbance.  As  a  result,  the  infinitesimal  three-dimensional  disturbance  was  rendered 
stable.  Clearly,  the  linear  controller  not  only  stabilized  the  two-dimensional  lin¬ 
ear  system,  but  also  had  a  stabilizing  effect  on  the  three-dimensional  nonlinear 
^vstem.  These  results  show  the  promise  of  linear  controllers  even  in  nonlinear  sys¬ 
tems. 


10.  Conclusions 

In  this  paper,  we  have  developed  feedback  controllers  that  linearly  stabilize 
plane  Poiseuille  flow.  We  used  a  Galerkin  spectral  method  to  generate  state-space 
models.  Indeed,  any  (convergent)  numerical  method  that  reduces  the  governing  par¬ 
tial  differential  equations  into  a  set  of  ordinary-differential  equations  may  be  used. 
It  is  important  to  note,  however,  that  the  meaning  of  the  state  variables  in  the  state- 
space  model  changes  as  different  numerical  methods  are  employed.  Even  though 
state  variables  may  not  have  any  specific  physical  meaning,  some  numerical  meth¬ 
ods  may  result  in  favourable  divisions  of  system  dynamics.  In  the  plane  Poiseuille 
flow  case,  the  spectral  Galerkin  method  with  the  use  of  Fourier  components  in  the 
x-direction  and  combination-Chebyshev  polynomials  in  the  y-direction  led  to  mod¬ 
elled  dynamics  that  were  de-coupled  by  wavenumber.  This  led  to  a  block  diagonal 
form  of  the  A  matrix.  As  a  result,  we  were  clearly  able  to  describe  modelled  and 
unmodelled  dynamics  in  terms  of  wavenumber  dynamics  included  and  not  included 
in  our  finite-dimensional  model.  This  also  led  us  to  the  concept  of  using  distributed 
control  to  render  certain  wavenumbers  ‘uncontrollable’.  If  other  numerical  methods 
had  been  used,  this  concept  would  not  have  been  so  transparent.  Furthermore, 
plane  Poiseuille  flow  can  be  linearly  stabilized  with  simple  feedback  controllers  if 
sensors  are  pla.ced  at  judicious  locations.  Both  the  position  and  the  type  of  sensing 
and  actuation  change  the  zeros  of  single-input/single-output  models.  In  our  case, 
we  have  seen  that  certain  shear  sensor  locations  lead  to  ‘minimum-phase’  systems 
and  some  locations  lead  to  ‘non-minimum  phase’  systems.  Minimum-phase  systems 
are  in  general  easier  to  control  than  non-minimum  phase  systems.  As  a  result,  we 
svere  able  to  stabilize  plane  Poiseuille  flow  with  a  constant  gain  feedback,  integral 
compensator  controller.  In  addition,  the  controller  was  extremely  robust  to  a  wide 
range  of  Reynolds  numbers.  Also,  we  have  shown  that  one  danger  of  feedback 
control  is  that  the  linear  transient  response  of  a  controlled  system  can  lead  to  high 
amplitudes  for  a  short  period  of  time.  If  these  amplitudes  are  high  enough,  it  is 
possible  that  they  may  invalidate  the  linear  model  and  enhance  nonlinear  effects. 
Furthermore,  the  high  transients  may  not  be  observable  at  the  output  so  that  feed¬ 
back  control  cannot  be  used  to  suppress  them.  Finally,  we  have  shown  that  linear 
controllers  have  a  strong  stabilizing  effect  on  two-dimensional  finite-amplitude  dis¬ 
turbances.  As  a  result,  three-dimensional  secondary  instabilities  can  be  rendered 
stable. 
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A  linear  process  in  wall-bounded  turbulent  shear  flows  has  been  investigated  through  numerical 
experiments.  It  is  shown  that  the  linear  coupling  term,  which  enhances  non-normality  of  the 
linearized  Navier-Stokes  system,  plays  an  important  role  in  fully  turbulent— and  hence,  nonlinear 
— flows.  Near-wall  turbulence  is  shown  to  decay  without  the  linear  coupling  term.  It  is  also  shown 
that  near-wall  turbulence  structures  are  not  formed  in  their  proper  scales  without  the  nonlinear  terms 
in  the  Navier- Stokes  equations,  thus  indicating  that  the  formation  of  the  commonly  observed 
near-wall  turbulence  structures  are  essentially  nonlinear,  but  the  maintenance  relies  on  the  linear 
process.  Other  implications  of  the  linear  process  are  also  discussed.  ©  2000  American  Institute  of 
Physics.  [S 1070-663 1  (00)00708-X] 


The  transient  growth  due  to  non-normality  of  the  eigen- 
modes  of  the  linearized  Navier-Stokes  (N-S)  equations  has 
received  much  attention  during  the  past  several  years  (see, 
for  example,  Refs.  1-3).  It  has  been  shown  that  the  energy  of 
certain  disturbances  can  grow  to  G{ Re2)  in  time  proportional 
to  O(Re),  where  Re  denotes  Reynolds  number  of  the  flow.2 
It  has  been  postulated  that  this  transient  growth,  which  is  a 
linear  process,  can  lead  to  transition  to  turbulence  at  a  Rey¬ 
nolds  number  smaller  than  the  critical  Reynolds  number,  be¬ 
low  which  a  classical  linear  stability  theory  based  on  the 
modal  analysis  predicts  that  all  small  disturbances  decay 
asymptotically.  As  such,  some  investigators  attributed  this 
linear  process  as  a  possible  cause  for  subcritical  transition  in 
some  wall-bounded  shear  flows,  such  as  plane  Poiseuille 
flow  and  Couette  flow. 

Some  investigators  further  postulated  that  the  same  lin¬ 
ear  process  is  also  responsible  for  the  observed  wall-layer 
streaky  structures  in  turbulent  boundary  layers.1,4  The  opti¬ 
mal  disturbance,  which  has  the  largest  transient  growth  ac¬ 
cording  to  their  optimal  perturbation  theory,  looks  similar  to 
the  near-wall  streamwise  vortices  that  create  the  streaky 
structures  in  turbulent  boundary  layers.  However,  this  opti¬ 
mal  disturbance  occupies  the  entire  boundary  layer,  in  con¬ 
trast  to  the  streamwise  vortices  in  turbulent  boundary  layers, 
which  are  confined  to  the  near-wall  region.  In  order  to  relate 
their  optimal  perturbation  theory  to  those  structures  observed 
in  turbulent  boundary  layers,  a  time  scale  corresponding  to 
the  bursting  process  in  turbulent  boundary  layers,  which  is 
essentially  a  nonlinear  process,  was  introduced  as  an  addi¬ 
tional  parameter.4  It  has  been  argued  that  the  transient 
growth  in  turbulent  boundary  layers  would  be  disrupted  by 
turbulent  motions  on  a  time  scale  corresponding  to  the  burst¬ 
ing  process,  which  is  smaller  than  the  viscous  time  scale,  and 


hence,  the  globally  optimal  disturbance  will  never  have  a 
chance  to  grow  to  its  maximum  possible  amplitude.  The  no¬ 
tion  that  commonly  observed  wall-layer  structures  are  related 
to  a  linear  process,  although  it  is  the  nonlinear  process  that 
determines  the  proper  length  scale,  suggests  that  the  same 
linear  process  may  play  an  important  role  in  fully  nonlinear 
turbulent  boundary  layers. 

Other  evidence  that  a  linear  process  may  play  an  impor- . 
tant  role  in  turbulent  boundary  layers  can  be  found  in  the 
work  of  Joshi  et  a/.5’6  and  others,7’8  who  successfully  applied 
controllers  developed  based  on  a  liner  system  theory  to  the 
nonlinear  flow  in  their  attempt  to  reduce  .he  viscous  drag  in 
turbulent  boundary  layers.  Bewley9  anr^ied  linear  optimal 
control  theory  to  a  nonlinear  convection  problem.  Although 
it  is  not  clear  how  controllers  based  on  a  linearized  model 
work  so  well  for  nonlinear  flows  and  it  is  a  subject  of  further 
investigation,  these  results  suggest  that  the  essential  dynam¬ 
ics  of  near-wall  turbulence  may  well  be  approximated  by  a 
linear  model. 

Motivated  by  the  above  findings,  we  investigate  the  role 
of  this  linear  process  in  fully  nonlinear  turbulent  flows.  In 
particular,  we  investigate  the  role  of  the  linear  coupling  term 
(see  below  for  its  definition),  which  is  a  source  of  the  non¬ 
normality  of  the  eigenmodes  of  the  linearized  Navier-Stokes 
equations,  in  wall-bounded  shear  flows,  using  a  turbulent 
channel  flow  as  an  example. 

In  this  Letter,  we  shall  use  (x,y,z)  for  the  streamwise, 
wall-normal,  and  spanwise  coordinates,  respectively,  and 
(u,v,w)  for  the  corresponding  velocity  components.  Rey- 
nolds  number,  ReT,  is  based  on  the  wall-shear  velocity,  uT 
-  Vt^/p,  and  the  channel  half-width,  h,  where  rw 
=  vdUldy\w  is  the  mean  shear  stress  at  the  wall,  and  U,  v , 
and  p  denote  the  mean  velocity,  viscosity,  and  density,  re- 


1 070-6631/2000/1 2(8)/1 885/4/SI  7.00 


1885 


©  2000  American  Institute  of  Physics 


1886  Phys.  Fluids,  Vol.  12,  No.  8,  August  2000 


J.  Kim  and  J.  Urn 


spectively.  The  superscript  “  +  ”  denotes  quantities  nondi- 
mensionalized  by  v  and  u  T . 

Representing  the  wall-normal  velocity,  v,  and  the  wall- 
normal  vorticity,  o)y)  in  terms  of  Fourier  modes  in  the 
streamwise  ( x )  and  the  spanwise  (z)  directions,  the  linear¬ 
ized  N-S  equations  can  be  written  in  an  operator  form 

d  v  v 

T.  ,  .  ,  (1) 


and  the  hat  denotes  a  Fourier-transformed  quantity.  Here 
Los,  £Sq>  and  Lc  represent  the  Orr-Sommerfeld,  Squire,  and 
the  coupling  operators,  respectively,  and  defined  as 

L05=A-](-ikxim  +  ikx(d2U/dy2)  +  (l/Re)A2), 

Isq=-i^C/+(l/Re)A,  (3) 

Lc=  —  ikz  (dUldy) , 


where  kx  and  k,  are  the  streamwise  and  spanwise  wave  num¬ 
bers,  respectively,  k2  =  k2x  +  k2  ,  A  =  d1ldy1—k2,  and  U  is  the 
mean  velocity  about  which  the  N-S  equations  are  linearized. 
Note  that  the  full  nonlinear  N-S  equations  can  be  written 
also  as 


d 

dt 


(4) 


where  all  nonlinear  terms  are  lumped  into  Mv  and  The 
operator  A  in  this  case,  however,  is  a  function  of  v  and  cuv , 
because  U  depends  on  u  and  o>v . 

It  has  been  shown  that  operator  A  in  Eq.  (2)  is  non¬ 
normal,  and  hence,  its  eigenmodes  are  nonorthogonal,  thus 
allowing  a  transient  growth  of  energy  even  if  all  individual 
modes  are  stable  and  decay  asymptotically.1,3  Note  that  the 
coupling  term  Lc  vanishes  for  two-dimensional  (2-D)  distur¬ 
bances  (£z  =  0),  and  therefore,  there  is  no  coupling  between 
v  and  tov  for  2-D  disturbances.  For  3-D  disturbances,  how¬ 
ever,  v  evolves  independently,  but  cov  is  forced  by  v  through 
the  coupling  term.  It  should  be  noted  that  Los  itself  is  not 
self-adjoint,  and  hence,  2-D  disturbances  can  have  a  transient 
growth,  but  it  was  shown  that  3-D  disturbances  have  much 
larger  transient  growth  due  to  the  coupling  term,  which 
causes  larger  non-normality.  In  the  present  study,  we  concen¬ 
trate  on  the  role  of  the  coupling  term  in  fully  nonlinear  tur¬ 
bulent  flows,  using  a  fully  developed  turbulent  channel  flow 
as  an  example. 

In  order  to  investigate  the  role  of  the  coupling  term  in 
fully  turbulent  flows,  we  proceed  to  solve  the  following 
modified  nonlinear  equations: 


d  v  L0s  0  v 

toy  .0  ^sq  -  toy 

This  modified  system  can  be  viewed  as  representing  a  syn¬ 
thetic  turbulent  flow  without  the  coupling  term,  or  a  turbu¬ 


t* 

FIG.  1.  Time  evolution  of  mean  shear  at  wall: - ,  upper  wall;  , 

lower  wall.  Thick  lines  are  for  a  regular  channel  flow,  while  thin  lines  are 
for  a  channel  flow  with  Lc=  0  in  the  upper  half  of  the  channel  starting  from 
r+  =  0. 

lent  flow  with  control  by  which  the  coupling  term  is  sup¬ 
pressed.  For  instance,  surface  blowing  and  suction  activated 
to  eliminate  (reduce)  the  spanwise  variation  of  v  (i.e., 
dvfdz)  could  eliminate  (reduce)  the  effect  of  the  coupling 
term. 

A  spectral  channel  code  similar  to  that  of  Kim  et  a/.10 
was  used  to  solve  the  above  modified  nonlinear  equations. 
To  further  contrast  the  role  of  the  coupling  term,  we  used  the 
modified  N-S  equations  only  in  the  upper  half  of  the  channel 
and  the  regular  N-S  equations  in  the  lower  half  of  the  chan¬ 
nel.  We  used  the  same  Reynolds  number  (Rer=100)  and 
grid  (32X  65X32  in  x,y,z)  as  Lee  et  aLu 

In  the  first  numerical  experiment,  we  used  a  regular  tur¬ 
bulent  velocity  field  obtained  by  Lee  et  al  as  our  initial  field. 
Starting  from  this  initial  field,  we  integrated  in  time  to  see 
how  the  turbulent  flow  in  the  upper  half  of  the  channel 
evolves  in  the  absence  of  the  coupling  term.  Time  evolution 
of  the  mean  shear  at  both  walls  is  shown  in  Fig.  1,  which 
illustrates  a  drastic  reduction  in  the  wall  shear  without  the 
coupling  term.  Several  snapshots  of  the  velocity  field  are 
shown  in  Fig.  2,  where  contours  of  streamwise  vorticity  in  a 


FIG.  2.  Contours  of  streamwise  vorticity  in  y-z  plane:  (a)  t+  =  0;  (b)  /+ 
-20;  (c)  r+  =  200.  -80<wjr<80  with  18  contour  levels.  Note  that  Z,c=0 
only  in  the  upper-half  of  the  channel. 
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=  180. 

y  —  z  plane  are  shown  to  illustrate  the  effect  of  the  coupling 
term  on  turbulence  structures.  It  is  evident  that  streamwise 
vortices  quickly  disappear  without  the  coupling  term.  The 
reduction  of  the  wall  shear  in  conjunction  with  the  disappear¬ 
ance  of  the  streamwise  vortices  is  a  common  feature  of  many 
drag-reduced  turbulent  flows.11  Turbulence  intensities  shown 
in  Fig.  3  indicate  drastic  reductions  without  the  coupling 
term. 

In  the  second  numerical  experiment,  we  used  an  initial 
velocity  field  consisting  of  the  same  mean  velocity  as  the 
first  experiment  but  with  random  disturbances,  and  hence, 
there  are  no  organized  turbulence  structures  present  initially. 
A  divergence-free  white-noise  spectrum  was  used  for  this 
purpose.  The  amplitude  was  chosen  such  that  neither  they 
decay  too  quickly  (too  small)  nor  they  cause  a  numerical 
instability  due  to  non-smoothness  of  the  initial  condition  (too 


FIG.  4.  Contours  of  streamwise  vorticity  in  y  —  z  plane  at  C  =  20,  started 
from  an  initial  random  field:  (a)  Case  1,  regular  turbulent  flow;  (b)  Case  2, 
without  the  linear  coupling  term,  iz\  (c)  Case  3,  without  the  nonlinear 
terms.  Contour  levels  are  the  same  as  Fig.  2. 


A  linear  process  in  wall-bounded  turbulent  shear  flows  1887 


FIG.  5.  Contours  of  streamwise  vorticity  in  y—z  plane  at  r+=40,  started 
from  a  random  initial  field.  See  figure  caption  in  Fig.  4  for  legend. 

large).  Starting  with  the  same  random  initial  field,  three  dif¬ 
ferent  simulations  were  carried  out:  Case  1,  with  the  full 
nonlinear  N-S  equations  (i.e.,  regular  turbulent  flow);  Case 
2,  with  N-S  equations  without  the  linear  coupling  term; 
Case  3,  with  N-S  equations  without  the  nonlinear  terms  (i.e., 
linearized  N-S).  The  purpose  of  these  simulations  is  to  in¬ 
vestigate  whether  the  linear  coupling  term  is  indeed  respon¬ 
sible  for  formation  of  the  streamwise  vortices  and  near-wall 
streaks,  and  if  so,  whether  the  time  scale  associated  with  the 
formation  of  these  structures  corresponds  to  the  bursting  pro¬ 
cess  (/  +  ~100),  as  hypothesized  by  Butler  and  Farrell.4 

Time  evolution  of  the  three  velocity  fields  is  shown  in 
Figs.  4-6  with  streamwise  vorticity  contours  in  a y-z  plane. 
Organized  structures  are  discernible  in  all  three  cases  as 
early  as  f+-20  (Fig.  4),  but  they  look  different  from  each 
other.  For  Case  3  (without  the  nonlinear  terms),  the  struc¬ 
tures  that  appear  from  the  structureless  random  initial  condi¬ 
tion  have  larger  spanwise  scales  than  those  in  the  regular 
flow  (Figs.  4  and  5).  For  Case  2  (without  the  linear  coupling 
term),  the  vortical  structures  appear  briefly  (Fig.  4)  but  dis¬ 
appear  quickly,  especially  in  the  wall  region  (Figs.  5  and  6), 
since  they  cannot  be  maintained  without  the  linear  coupling 
term  as  demonstrated  in  the  first  experiment  mentioned 
above.  For  Case  1,  the  time  for  these  structures  to  appear  is 
shorter  than  that  implied  by  the  optimal  disturbance  mecha¬ 
nism  of  Butler  and  Farrell.4  Note  that  the  structures  in  Case  3 
are  already  substantially  different  from  those  in  Case  1  at 
r+  =  40,  indicating  that  the  effect  of  nonlinear  terms  is  felt 
much  earlier  than  the  eddy  turnover  time  (t+^  100)  as  pro¬ 
posed  by  Butler  and  Farrell.4  The  present  result  is  also  con¬ 
sistent  with  Jimenez  and  Pinelli,12  who  showed  that  the  for¬ 
mation  of  streaky  structures  can  be  prevented  by  damping 
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FIG.  6.  Streamwise  contours  in  y  —  z  plane  at  r+  =  80,  started  from  a  ran¬ 
dom  initial  field.  See  figure  caption  in  Fig.  4  for  legend. 

(vo):)  [the  (  )  indicates  streamwise  average],  which  is  re¬ 
lated  to  the  linear  coupling  term.  While  both  results  clearly 
demonstrate  the  essential  role  of  the  linear  coupling  term  in 
the  formation  and  maintenance  of  the  wall-layer  streaks, 
while  the  present  work  also  indicates  that  a  nonlinear  mecha¬ 
nism  is  responsible  for  producing  the  proper  streak  spacing. 

We  have  used  several  different  initial  conditions  to  de¬ 
termine  whether  the  above  results  depend  on  initial  condi¬ 
tions,  but  found  no  such  evidence.  It  thus  appears  that  both 
the  nonlinear ierms  and  the  linear  xouplingjvim  are  neces¬ 
sary  for  the  formation  and  maintaining  of  these  structures  at 
their  proper  scale.  The  nonlinear  terms  are  necessary  for  the 
formation  of  streamwise  vortices  and  the  linear  coupling 
term  is  necessary  to  generate  the  wall-layer  streaks,  the  in¬ 
stability  of  which  in  turn  strengthen  the  streamwise  vortices 
through  a  nonlinear  process.  In  the  absence  of  either  mecha¬ 
nism,  turbulence  ceases  to  exist.  The  result  of  this  second 
experiment  is  consistent  with  Hamilton  et  aln  and  Waleffe 
and  Kim14  in  that  the  formation  of  the  streamwise  vortices  is 
a  result  of  a  nonlinear  process. 

We  have  demonstrated  that  the  linear  process  associated 
with  the  coupling  term  plays  an  important  role  even  in  fully 
nonlinear  wall-bounded  turbulent  shear  flows.  Near- wall 
streamwise  vortices,  which  play  the  essential  role  in  the  dy¬ 
namics  of  wall-bounded  shear  flows,  are  seen  to  be  formed 
but  cannot  be  sustained  without  the  coupling  term. 

This  result  is  consistent  with  the  analysis  by  Henningson 
and  Reddy15  who  showed  that  non-normality  of  the  linear¬ 
ized  Navier-Stokes  operator  is  a  necessary  condition  for  dis¬ 
turbances  to  grow  for  Reynolds  number  below  the  critical 
Reynolds  number  predicted  by  the  traditional  linear  stability 
analysis.  However,  we  believe  this  is  the  first  direct  demon¬ 


stration  that  turbulence  (nonlinear  disturbance)  decays  when 
the  non-normality  of  the  underlying  linear  operator  in  non¬ 
linear  flows  is  reduced.  The  time  scale  associated  with  the 
formation  is  found  to  be  smaller  than  the  bursting  process 
used  in  the  optimal  perturbation  theory.  The  fact  that  the 
coupling  term  plays  an  essential  role  in  maintaining  the 
streamwise  vortices,  which  have  been  found  to  be  respon¬ 
sible  for  high  skin-friction  drag  in  turbulent  boundary  layers, 
suggests  that  an  effective  control  algorithm  for  drag  reduc¬ 
tion  should  be  aimed  at  reducing  the  effect  of  the  coupling 
term  in  the  wall  region.  In  fact,  the  opposition  control  used 
by  Choi  et  al 16  can  be  viewed  as  a  control  scheme  trying  to 
reduce  the  effect  of  the  coupling  term  by  suppressing  the 
span  wise  variation  of  v  in  the  wall  region.  It  should  be  in¬ 
teresting  to  design  a  control  algorithm  that  directly  accounts 
for  the  coupling  term  in  a  cost  function  to  be  minimized. 

A  portion  of  this  work  was  carried  out  while  the  first 
author  was  visiting  the  Isaac  Newton  Institute  for  Math¬ 
ematical  Sciences  at  Cambridge  University.  This  work  has 
been  supported  by  the  Air  Force  Office  of  Scientific  Re¬ 
search  (F49620-97- 10276,  Dr.  Marc  Jacobs).  The  computer 
time  has  been  provided  by  the  San  Diego  Supercomputer 
Center. 
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